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ABSTRACT 


Until recently, computer simulations of helicopter rotor dynamics have employed 
equations of motion that have been linearized or simplified. These modified equations of 
motion did not allow for the evaluation of nonlinear material properties in the rotor since 
higher order terms in the dynamics had been modified in the simplification process. With 
recent advances in both computer simulation hardware and symbolic mathematic 
manipulation software, the full nonlinear equations of motion may be utilized in helicopter 
rotor simulations. This dissertation reports on the use of the full nonlinear equations of 
motion in the analysis of rotor blade lead/lag motion and its effect on rotor hub and rigid 
body fuselage motion. Nonlinear modeling methods are implemented using Maple 
symbolic mathematic manipulation software and Matlab and Simulink computer 
simulation environments. Results are compared to the RAH-66 Comanche Froude scale 
wind tunnel article and new methodologies evaluated in the search for a damperless rotor 
system that is free of ground and air resonance mechanical instabilities. 
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I. INTRODUCTION 
One of the goals of this research was to develop a 
flexible computational tool to analyze the dynamic and 
aeromechanical behavior of advanced technology coupled 
rotor/fuselage systems. Initially, a series of programs were 


developed utilizing the symbolic processing software, 


Maple?, the computational software, Matlab”, and the 


simulation software, Simulink®, by LT Christopher Robinson 


[Ref. 1]. It was desired that the computational tool be 
simple to understand and lend itself to easy reprogramming 
by any user knowledgeable in the field of dynamics and in 
the use of the software programs mentioned. It was also 
desired that the developed programs allow for a sufficient 
capability so that the effects of introducing advanced 
technologies into rotor system designs could be accurately 
modeled. This dissertation reports on further advances in 
this computational tool development, and analysis based on 
these new tools. The nonlinear rotor simulation that was 
developed as part of this research is a very powerful tool 
for analyzing rotors that do not have multiple load paths to 
the hub. 

Historically, rotor analysis was performed using 
approximate equations of motion. In order to reduce the 


complexity of the computer code involved and to reduce 


computation times, these equations were simplified by 
eliminating higher order terms [Refs. 2,3]. Subsequently, 
either linearized equations were used or select nonlinear 
terms were retained using a ranking system or ordering 
scheme, such as that used by Friedmann, [Ref. 4]. 

As  hingeless helicopter main  rotors became more 
commonplace, a need arose for a rotor simulation tool that 
would accurately model nonlinear mechanical properties so 
that these nonlinearities could be exploited. It was desired 
to develop an analysis tool that would be able to reliably 
model the effects of nonlinearities in the rotor and hub 
mechanical and geometric parameters. For this reason, 
Robinson and Wood developed a utility for creating the full 


nonlinear, coupled equations of motion. Robinson utilized 
the symbolic manipulation software Maple? for this purpose. 


Current trends in helicopter technology and 
manufacturing have favored the use of bearingless rotor 
designs that make use of advanced composite materials. These 
designs offer many advantages over more conventional 
articulated rotors in reliability and maintainability. A 
potential payoff from the successful use of the technologies 
mentioned above is the damperless rotor; a design that 
offers major returns in the form of decreased rotor system 
weight, reduced parts count, and reduced maintenance 


requirements. 


Composites and other advanced materials that can be 
applied for a damperless rotor make modeling rotor behavior 
more difficult, however. These materials have the potential 
for exhibiting nonlinear behavior that cannot be accurately 
modeled without utilizing the full nonlinear equations of 
motion for the system. 

The goal of this dissertation is to examine alternative 
designs for a damperless helicopter rotor. This effort met 
with success by taking advantage of nonlinear stiffness in 
the blade root end to avoid divergent motion in the rotor 
blade lead/lag degrees of freedom. This oscillatory 
instability is better known as ground and air resonance 


instability. 








II. REVIEW OF WORK PERFORMED BY ROBINSON 
For background on progress reported iji this 
dissertation, it is important to first review LT Robinson's 


thesis work. This work preceded that of the author. Tasks 
performed by Robinson included the formulation of a Maple? 


based symbolic processing worksheet that formulated 


nonlinear eguations of motion given energy expressions for 
helicopter rotor model degrees of freedom. Simulink? based 


computer simulations were developed from the equations of 
motion derived by the symbolic processor for a simple three 
bladed rotor based on that used by Coleman [Refs. 5,6]. 
Robinson, Wood, and King reported on a new method for 
formulating the full non-linear equations of motion for 
ground/air resonance stability analysis of helicopter rotor 
systems [Ref. 7]. The full set of non-linear equations was 
developed by Lagrangian approach using symbolic processing 
software for expanding the equations. The symbolic software 
was further utilized to automatically convert the equations 


of motion into C or Fortran source code formatted for 
numerical integration. Simulink? then applied a Runga-Kutta 


integration scheme to generate time history plots of blade 


and fuselage motion. Damping levels were determined from the 


time history simulations by a Matlab? program, which used 


the Moving Block Technique for establishing damping levels 
from the coupled rotor-body response. 
Robinson's model, based on that of Coleman, is shown in 


Figure 1, [Ref- T]. 





Figure 1 - Simplified Rotor Model [After Ref. 1] 


Robinson verified his symbolic Lagrangian derivation by 
comparing his results to Coleman's. Figure 2 shows results 
of a comparison of the two methods for an unstable case 
where a moderate amount of damping has been added to both 


rotor blades and fuselage. 


Comparison of Simulation Results to Solution of Coleman Equations 


. [a] 
Simulation 
s Coleman 


o 
u 


2 
S 
O 
o 
= 
TE 
o 
E 
o 
O 
& 
o. 
2 
O 
Oo 
9 
O 
o 
o 
SH 


4 
— 


time (sec) 





Figure 2 - Validation of simple model with solution of Coleman's linearized model 


In Figure 2, Robinson shows excellent agreement between 
the two solutions with departure occurring only when 
displacements are very large. This is to be expected since 
the simulation model of this paper does not assume small 
angle theory whereas the Coleman-Feingold-Bramwell model 
does. 

For a more quantitative proof of the validity of the 
new method, it was tested against Deutsch's Criteria (Ref. 
8]. Based on Coleman's analysis, Deutsch showed that the 
product of the hub (landing gear) damping and the blade 


lead/lag damping must be greater than a prescribed value to 


collapse the band of unstable blade frequencies: A, Ag > A 


(p-1). The following Figure shows results of successive 
sweeping across the unstable band using the new nonlinear 
analysis. The respective cases represent those where 
Deutsch’s prescribed damping is applied (D=1.0), and other 
cases where the damping is greater than that required by 
Deutsch (D=1.1) or less (D=0.8, 0.9). For each case 
considered, the center of instability is seen to be the 
point of minimum value on the given line. With the new 
nonlinear analysis, Deutsch’s criteria applied at the center 
of instability provides nearly a neutrally stable result 
with a near zero damping ratio. Deutsch’s criteria is 
conservative in that the critical point for his criteria, or 
bucket of each curve, still shows positive damping when 
analyzed with the NPS simulation. Three years after 
Robinson, Rafanello verified the conservative nature of 
Deutsch’s criteria by matching simulation results to Wood’s 


HSS-2 ground resonance analysis, [Refs. 9, 10]. 
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Figure 3 - Quantitative validation of simple model against Deutsch Criteria 








III. DAMPERLESS ROTOR DEVELOPMENT 


A. BACKGROUND 

Ground and air resonance are particularly destructive 
mechanical instabilities that can occur in helicopters with 
fully articulated, bearingless, or hingeless main rotor 
designs. The phenomenon of ground/air resonance is the 
result of coupling between rigid body motions of the 
fuselage on its landing gear, or in the air, and lead-lag 
oscillations of the rotor blades in their plane of rotation. 
When it occurs, the instability can build to destructive 
proportions in a matter of seconds. 

The equations of motion describing the dynamics of the 
coupled rotor-fuselage system are nonlinear and generally 
quite complex even for simplified models. The fundamental 
theory of air and ground resonance is credited to the 
classic work of an NACA engineer, Robert Coleman, who first 
published his theory in the early 1940's (Ref. 5]. As 
computational power improved with the evolution of digital 
computers, more general techniques for analyzing rotor 
system stability were developed. 

As stated earlier, due to the complexity of the 
equations of motion, ground and air resonance analysis was 
historically restricted to linear analysis until very 


recently. This frequently restricted the scope of many 


PL 


studies and potential solutions. With a few noteworthy 
exceptions, nonlinearities in ground and air resonance 
research programs have not typically investigated beyond 
traditional damping methods, such as hydraulic damping. 

This chapter entitled "Damperless Rotor Development" 
investigates various options for damperless rotors that 
maintain some margin of stability, or avoidance of ground 
and air resonance without the employment of auxiliary 
lead/lag dampers. 

In general, rotorcraft lag dampers are heavy, expensive 
items that can require extensive maintenance. By developing 
rotors that do not require these dampers, payoffs in 1) 
reduced maintenance time and cost, 2) reduced hub complexity 
and parts count, 3) reduced rotor drag, and 4) reduced rotor 
signature could potentially be reaped. 

An effort in the elimination of these dampers applies 
directly to a famous quote from the father of the business 
jet, Bill Lear, "Strive for design simplicity: you never 


have to fix anything you leave out." 


B. GROUND RESONANCE FUNDAMENTALS 

Ground and air resonance instability is centered around 
blade/hub interactions at one of the fuselage rigid body 
nodes. The "rotor" acts as a forcing function, producing 


undesired responses or even resonance in the fuselage 


12 


dynamics. The rotor changes the frequency of mechanical 
vibrations as they transfer from the rotating frame down 


into the non-rotating frame. If the blades have a primary 


mode in lead/lag at some frequency 0 and 2 is the 


lag! 


frequency of rotation of the rotor system, shears are 
transmitted through the mast and into the non-rotating hub 
is 


at frequencies (Qto 


ine) B 


The progressing mode (Q+@ 


SEN 
heavily damped and normally plays no role in ground 


resonance. The  regressing mode  (Q-0,.) 


provides the 
potential for destruction, however. If one of the fuselage 
rigid body modes is near in frequency to the regressing mode 
and there is insufficient damping, blade lead/lag and 
fuselage displacement amplitudes can build to destructive 
proportions. While the helicopter is still in contact with 
the ground, the fuselage roll mode is often close in 
frequency to the regressing lag mode, normally in the 2-5 Hz 
range. 

In general, the spectrum of classical ground resonance 
to air resonance could be referred to as "aeromechanical 
stability.” Since articulated rotors “normally do not 
experience air resonance, aerodynamics and flap degrees of 
freedom play little roll in ground resonance. Therefore, 
articulated rotors may be modeled using Coleman's methods. 
Hingeless and bearingless rotors experience more coupling in 


the blade degrees of freedom. In addition, body pitch and 


T> 


roll coupling through flap moments of cantilever blades can 
result in air resonance. Bearingless and hingeless rotors, 
where Coleman's equations are not sufficient to model the 
problem, require more sophisticated models. Flap degrees of 
freedom and aerodynamics were not used to model the Coleman- 
like rotor models in this report. 

In the typical rotor, there are, in fact, significant 
nonlinearities. Especially noteworthy are friction in 
dampers and oleos that sometimes result in the appearance of 
instability following a moderately large external force 
excitation. Thus, the typical scenario for ground resonance 
is initiated by the fuselage being perturbed by some 
external force such as a gust, a rolling ship deck, uneven 
terrain while ground taxing, or combat damage. The fuselage 
motion, typically in roll, induces asymmetrical lag motion 
in the rotor blades. This asymmetrical response shifts the 
combined CG of the rotating components away from the center 
of rotation, resulting in an unbalance in the centrifugal 
force seen by the rotor mast. The unbalanced centrifugal 
force results in increased fuselage motion, compounding the 
problem. This cycle is depicted in flowchart form in Figure 


4. 
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Figure 4 - Ground resonance flowchart 


A typical fuselage on its landing gear has a roll mode 
frequency in the proximity of the regressing mode frequency 
of the rotor system. If these modes are sufficiently near to 
each other and there is insufficient system damping, ground 


resonaneec= cal Cecut. 


m5 


A plot of hub motion for a typical case of ground 
resonance can be found in Figure 5. Note the spiraling 
divergence of the hub center of mass in the horizontal 


plane. 


Hub Motion in the Horizontal Plane 
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Figure 5 - Ground resonance hub motion 


Following the work of Robert Coleman, a more complete 
work, based on his analysis, was collected and published by 
Coleman and Feingold [Ref. oj: Deutsch supplemented 
Coleman’s work by establishing a criterion for system 


damping necessary to eliminate ground resonance [Ref. 8]. 
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Numerous other researchers have contributed EO 
understanding the ground resonance problem since Coleman, 
Deutsch, and Feingold. The papers following Coleman 
essentially completed the investigation of the linearized 
system for ground resonance. Hammond applied Floquet theory, 
but to consider the case of unequal lead/lag damping among 
the blades [Ref. 11]. (The NPS rotor simulation can also 
apply unequal lead/lag damping values at each blade.) 
Tongue, Flowers, Jankowski, Tang, Dowell, and most recently, 
Gandhi and Chopra  [Refs. 12-14] investigated nonlinear 
effects in blade and fuselage damping. Ormiston explored 
linear aeromechanical stability of rigid blades with spring 
restrained hinges, and later, investigated rotor modeling 
with a sophisticated finite element analysis code and 


various nonlinear damping models [Refs. 15, 16]. 


c. OBJECTIVES 

The objective of the present work is to explore the 
potential of eliminating the snubber-damper or damper on 
hingeless rotor designs and replace it with a flexbeam that 
has been modified to possess nonlinear properties. Thus, the 
purpose is to employ the nonlinear properties of the 
flexbeam to introduce nonlinear dynamic characteristics that 
will replace unbounded instability with small amplitude 


mit cycle oscillations, and prevent catastrophic 
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destruction of the helicopter. This work is well suited to 
the new NPS rotor simulation analysis that can accurately 
model nonlinear mechanical properties. 

The significance of this new research is that allowance 
of nonlinearities at the blade root may result in an 
acceptable bounded response in the parameter region where 
linear theory would predict instability, or may introduce 
new regions of unacceptable response in the parameter region 
where linear theory would predict stability. Evidence of the 
latter may be found in the numerous aircraft lost in the 
documented cases of ground resonance. Modern soft-inplane 
rotors, such as that first introduced on the BO-105, have 
the additional possibility of encountering this lead lag 
INstap Mita eh "REF. 17]. 

The work reported here attempts to abandon the use of 
auxiliary damping methods in the rotating system as a means 
of eliminating ground and air resonance. Only the damping 
inherent to the typical landing gear system and the built-in 
structural damping of the flexbeam in bending are utilized 
in the cases shown. The ground/air resonance cycle is 
detuned by shifting the blade natural frequency in lead/lag 
to avoid the coalescence of the regressing mode with the 
fuselage rigid body modes. A nonlinear  flexbeam is 
characterized by change in stiffness properties with change 


in amplitude. Therefore, as blade  lead/lag amplitude 
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changes, there is a corresponding change in its lead/lag 
frequency. This method differs from some previous methods of 
preventing ground resonance in that coalescence is avoided 
through natural frequency shift instead of dominating the 
dynamics of the system through damping. In the case of this 
thesis, this is accomplished by introducing a classic 
nonlinear stiffness known as Duffing-type nonlinear 


stiffness at the blade root. 


D. INTRODUCTION TO DUFFING 

L: Single Degree of Freedom Duffing Systems 

For the first part of the investigation, an exploration 
of the one degree of freedom Duffing system is considered to 
understand the effects of cubic stiffness for a simple 
harmonic oscillator. With damping excluded, Bufring s 


eguation can be expressed as: 
+: 2 Sn 
Xx + ox +hx~ = Gcos(0,!) 


pere the driving function is Gcos (0t): The coefficient iG, 
is the ratio between the amplitude of the force and the mass 
of the moving system. Coefficients o and h depend on the 


parameters of the system, with the Duffing coefficient, h, 
being positive for a "hardening spring" and negative for a 


Bsottening spring.” 
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One of the interesting phenomena appearing in systems 
governed by Duffing's equation is that, in addition to 
generating harmonics of the driving frequency, discontinuous 
jumps in amplitude occur as the frequency of the exciting 
force is either increased of decreased at a constant rate. 

A simple one degree of freedom Duffing spring mass 


system is shown in Figure 6. 


F * sin(w*t) 





Figure 6 - Simple Duffing system 


A simple simulation was written to obtain time 
histories for this system with user supplied mass, linear 


damping, linear stiffness, forcing amplitude, forcing 
frequency, and Duffing stiffness terms. The Simulink® 5” 


order Runge-Kutta integration scheme is used here. 
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An example time history that illustrates the effect of 
exciting frequency (both increase and decrease) is presented 


in Figure 7. 


Duffing System Time Histories 


Displacement (in) 
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Figure 7 - Duffing time history 


The inputs for this time history are contained in the 


following table: 


Table 1 - Table caption 
Linear stiffness = 4000 ib/in Mass = 6 slugs 


Duffing stiffness = 100 Ib/in Damping = 14 355.715 





ad 


Driving Force 1/2 amplitude = 500 lb 


The top subplot in Figure 7 as the time history of the 
mass displacement with a stiffening Duffing spring included 
in the dynamics. The second subplot in this figure is the 
driving frequency used in the top plot as a function of 
time. The system goes through resonance both as the driving 
frequency is increased and again as it decreases. 

If one takes this time history and transforms it as a 
function of driving frequency, instead of time, 
dissimilarities in the two resonances can easily be 


discerned, as shown in Figure 8. 
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Mass Deflection Time Histories 


Displacement (in) 





0 1 2 3 4 5 6 VA 8 9 10 
Frequency (Hz) 


Figure 8 - Duffing system resonance 


The system reaches resonance as the frequency is 
increased, but unlike a linear system that gradually reduces 
the displacement amplitude, in this case, the Duffing system 
experiences a sudden drop-off in displacement amplitude. 
Then, as the driving frequency is decreased from a value 
above resonance, the system does not follow the same path in 
resonance amplitude, but reaches a relative maximum in 


displacement at a frequency less than the case for 
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increasing driving frequency. These sudden changes in 
displacement amplitude are referred to as "Duffing jumps." 
Shown in Figure 9 is a nonlinear resonance curve for a 
hardening Duffing spring system, reproduced from Cunningham 
[Ref. 18]. A typical Linear resonance curve has been 
superimposed on Cunningham's original figure. The resonant 


frequency for the linear system would not change, but remain 
fixed at 0,. The physical reason for these jumps can be 


visualized in this plot of the resonance of the nonlinear 
Duffing system. 

As the driving frequency is increased from the left on 
the plot, the system displacement amplitude grows until the 
system reaches point A. Further increase in the driving 
frequency drops the displacement amplitude down to B, a 
Duffing jump. 

As the driving frequency is decreased from the right, 
the system drives right through point B to point C, where 
the solution becomes unstable. As the driving frequency is 
decreased lower than that found at point C, the only stable 
displacement amplitude becomes point D, another Duffing 


jump. 
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Figure 9 - Duffing resonance freguencies [After Ref. 18] 


Unlike the linear system, the Duffing system increases 
the natural frequency as a function of displacement 
amplitude. For larger displacements, the resonance frequency 
increases due to this stiffening. The relation for this 


frequency shift is derived by Cunningham as: 


[K 3K , A? 
0, <J—t —— 
M 4 


It should be noted that if the sign of h is allowed to 
reverse, h<1, (softening spring), the same kind of phenomena 


exists. The primary difference is that now the nonlinear 
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resonance curve is skewed towards the lower, rather than 
higher frequencies. Since ground and air resonance are 
determined by coupling of the rotor blades with lead/lag 
motion in the plane of rotation, it was hoped to take 
advantage of these sudden decreases in displacement 
amplitude. This decrease in blade lead/lag motion would 
subsequently reduce the severity of the offset centrifugal 
force on the hub. The employment of Duffing in the rotor 
lead/lag dynamics did, in fact, reduce the susceptability of 
the rotor to ground and air resonance, but not due to the 
Duffing jumps as hypothesized. The shift in natural 
frequency detuned the ground resonance sequence of events, 


instead. 


2% Multiple Degree of Freedom Duffing Systems 

Since the rotor system to be simulated was a three 
bladed design, investigations in a three degree of freedom 
Duffing system were also performed. A depiction of this 
system may be found in Figure 10. The possibility existed in 
this multiple DOF nonlinear system that reductions in 
resonant displacement amplitude could be achieved through 
coupling between the three degrees of freedom in the 


helicopter rotating frame. 
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be 


F * sin(w*t) 





Figure 10 - Three DOF Duffing system 


As shown, the Duffing spring/mass system described 
earlier was expanded into a three mass system. Each mass was 
linked by adjacent linear springs, linear dampers, and 
Duffing springs. An example time history of this three DOF 
system follows. The input parameters are the same as the one 
DOF case shown previously. The three resonant peaks are 


clearly visible. 
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Duffing System Time Histories 


15 20 25 
Time (sec) 


8 10 
Freq (Hz) 





Figure 11 - Three DOF Duffing system time histories 


In general, no substantial changes in the dynamics in 
the neighborhood of the first resonant frequency of this 
System were noted over the one DOF system. The second and 
third modes were most effected by variation in the Duffing 
terms of the system. By increasing the Duffing terms, the 


higher order modes tended to flatten out, showing resonant 
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behavior over larger ranges of the driving frequency. Three 
plots showing this trend are shown in Figure 12, first the 
linear system with no Duffing stiffness, with an 
intermediate Duffing coefficient, then with a large Duffing 


coefficient. 


Increasing Duffing Stiffness in 3 DOF Duffing System 
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Figure 12 - Effect of increasing Duffing stiffness in three DOF system 
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Since the ground and air resonance problem is more 
closely associated with the lower resonant frequencies of 
the rotor, no advantage in the ground/air resonance problem 
was gleaned from these higher order Duffing system analyses. 
The shift in resonance frequencies, however, still showed 
promise in rotor applications, as discussed in the following 


section. 


E. DUFFING ROTOR SYSTEMS 

As discussed previously, ground and air resonance are 
effectively driven by two phenomena: 1) the proximity of 
rotor and hub vibratory modes, and 2) the ability of system 
damping to maintain reasonable  DOF amplitudes when 
approaching resonance. Damping, the traditional means of 
solving the ground resonance problem, can dominate the 
system dynamics through the velocity terms. Except for high 
damping ratios, this approach does little to shift system 
modal frequencies, however. 

The method of passive rotor stability presented here 
utilizes a nonlinear Duffing type spring to address the 
coalescence part of the problem, while attempting to settle 
for whatever damping is inherent to the system. This 
approach eliminates the need for auxiliary damping devices 


for the blade lag degrees of freedom. 


30 


Coleman Model 2 Springs: 
| Duffing 
Linear 


Ng 
Structural 
Damping O 


Figure 13 - Rotor with Duffing stiffness in the lead/lag dynamics 





Duffing stiffness is added into the Coleman model as 
depicted in Figure 13. Both linear and nonlinear rotary 
springs plus structural damping are incorporated across the 
blade lag hinge. Stiffness and damping values can be defined 
separately for each blade to allow for asymmetric blade 
property investigations. 

The addition of a hardening Duffing spring increases 
the lead/lag natural frequency of the blade as a function of 
blade lead angle, detuning the resonance. As blade 
amplitudes increase in a typical ground resonance sequence 
of events, the regressing mode is driven to a lower 
frequency as the blade natural frequency in lag is increased 
due to stiffening. This effectively breaks the sequence of 


events previously depicted in flowchart form. This shift in 
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the rotor's regressing mode puts the blade lag motion in a 
limit cycle instead of letting displacements increase 
without bound. For the case of soft in-plane rotors that 
drive through the ground resonance region of instability 
upon rotor engagement, upon leaving this region, the 
dynamics return to what one would expect from a soft in- 
plane rotor at rotor speeds above the region of instability. 
Outside the region of instability, the linear terms once 
again dominate the dynamics, with little or no ill effects 


of the Duffing stiffness seen in the dynamics. 


1. Results of the Baseline Duffing Rotor 

The NPS nonlinear rotor simulation was utilized in the 
following investigation. To provide the most robust 
application of Duffing springs to the ground resonance 
problem, simulations are performed at the center of the 
region of instability unless otherwise stated. This baseline 


case is of a rotor with the following values: 


Table 2 - Table caption 


Q = 25.0 rad/s M,, < 180.0 sl Ma. = 6.0 sl 
(239 RPM) (5800 lbs) (193 lbs) | 


Structural Damping in Lag = 2.0% 
Viscous Damping in Hub = 2.0% 
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Time histories of blade and hub displacements for this 


rotor with no Duffing stiffness employed are illustrated in 


Figure 14 and Figure 15. 


The first plot shows the blade lag 


diverging 


initlal zero 


their 


from 


clearly 


angles 


displacement angles. 


Baseline Case w/o Duffing Spring 
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Figure 14 - 


The X direction of the isotropic hub motion is plotted 


in Figure 15. The hub has an initial lateral displacement of 


0.2 feet to initiate the ground resonance behavior of the 


Bo 


system. Like the blade plot above, the hub motion clearly is 


subject to an oscillatory divergence. 


Baseline Case w/o Duffing Spring 





Hub Displ. (ft) 


Time (sec) 


Figure 15 - Baseline case, hub lateral displacement without Duffing stiffness 


2. Baseline Case with Duffing Stiffness Added 

The next set of time histories gives results for the 
same rotor, but with Duffing spring terms added. The ratio 
of the Duffing stiffness coefficient to the linear stiffness 


coefficients 20.0 for Figure 16 and Figure 17. 
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Baseline Case w/ Duffing Spring 
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- Baseline case, blade lead angles with Duffing stiffness 


Figure 16 


The nonlinear stiffness clearly keeps the previously 


Since the rotor for 


unstable blade lag response in check. 


cycle of 


Tamie 


a 


unstable, 


case is linearly 


this 


approximately 8 degrees in the blade lag motion is produced. 


Limit cycle motion, 


A plot of the hub motion follows. 


is not considered desirable because it results 


in general, 


structural vibration 


additional 


non-rotating 


the 


in 


in 


The amplitude of the lateral motion 


aircraft components. 
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depicted in this plot is sufficiently large that concerns 


over hub vibration levels were raised. 


Baseline Case w/ Duffing Spring 
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Figure 17 - Baseline case, hub lateral displacement with Duffing stiffness 


The lateral vibrations experienced by the hub are 
plotted in Figure 18. The derivative of the hub velocity was 
taken numerically to obtain the hub lateral vibration data 
plotted. Since numerical derivation of a sampled signal 
tends to increase the effect of noise in the data, this plot 


is provided as a means of discerning the overall vibration 
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levels, not as an explicit means of calculating the 
vibration time history. The final vibratory levels are on 
the order of one-half the acceleration of gravity. 
Vibrations of this magnitude might be larger than desired, 


even for transient ground operations. 


Baseline Case w/ Duffing Spring 
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Figure 18 - Hub lateral vibration, baseline case with Duffing stiffness 


ST 


3% Baseline Case with Increased Rotor Speed - Above 
the Region of Instability 


The case that produced the unacceptable vibration 
levels presented above is one that was originally unstable 
and represents a worst-case operating regime for the rotor: 
the center of the region of instability. Even if the 
nonlinear stiffness was included in the rotor design, it is 
doubtful a production rotor would be designed to operate at 
its center of instability. 

Most soft in-plane rotors are designed to operate with 
rotor speeds slightly above the region of instability. The 
rotor drives through the unstable region on engagement and 
disengagement. It is also possible to encounter these low (2 
regimes during extreme maneuver, maintenance checks or as 
the result of a system failure. 

For the baseline rotor presented here, the blade center 
of mass is at a radial location of 10.5 feet; it roughly 
represents a 20-ft radius rotor. For this size blade, 25 
rad/sec (239 RPM) is an unusually low operating speed. If Q 
is increased to approximate a 650 ft/sec tip speed in a 
hover, (2 would increase to 32.5 rad/sec (310 RPM). 

To examine the effects of tHe inclusion of nonlinear 
stiffness at the higher rotor speed, plots of the blade and 
hub motion are provided in Figure 19 and Figure 20 for 


comparison with the linear rotor shown previously. 
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Base Case w/ Duffing & Incr Omega 
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Figure 19 - Blade lead angles, increased rotor speed with Duffing stiffness 


Since the rotor is now being operated with a rotor 


the hub and 


speed slightly above the region of instability, 


sama tial 


less For the same 


displacements are 


blade 


The lag limit cycle amplitude is reduced from 8 


eenditions. 


degrees all the way to 0.3 degrees. 
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Base Case w/ Duffing & Incr Omega 





Time (sec) 


Figure 20 - Hub lateral displacement, increased rotor speed with Duffing stiffness 


The hub lateral displacement saw a substantial 
decrease, as well. The hub limit cycle amplitude for the 
Duffing rotor drops 80$ from the value it had at the lower 
rotor speed. 

For comparison, the hub lateral vibrations for the 
Duffing rotor at both rotor speeds are plotted in Figure 21. 


The vibration amplitude drops almost an order of magnitude 
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a level that is certainly in the 


to less than 0.05 g"s, 


window of today's production rotorcraft. 


Base Case w/ Duffing & Incr Omega 
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Hub lateral vibration comparison 


Figure 21 - 


GROUND RESONANCE ANALYSIS 


F. 


Frequency Coalescence 


1. 


In the use of nonlinear stiffness to avoid ground and 


fundamental questions arise: 


In what way 


) 


d 


air resonance, 


does the stiffness effect the coalescence of the regressing 
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lag mode to the rigid body mode? and 2) Are the restoring 
moments produced of reasonable magnitude? 

Figure 22 is the Coleman Plot of the baseline rotor 
with and without Duffing stiffness included at the blade 
root end. It can be seen from the plot that the linear 
spring only case has coalescence with the hub lateral mode 
at a rotor operating RPM of 239, as depicted earlier in the 
divergent time histories, Figure 14 and Figure 15. With the 
inclusion of the stiffening Duffing spring, the non-rotating 
lag natural frequency is increased due to the increased lag 
displacement. This blade stiffening shifts the regressing 
lag mode away from the hub modal frequency at any operating 
RPM where resonance may occur. The action of the hardening 
Duffing spring keeps the frequency coelesence above the 


instantaneous rotor RPM. 
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Fixed System Fan Plot With Duffing Stiffness 
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Figure 22 - Coleman plot of baseline rotor with Duffing stiffness 


It should be noted that Duffing springs of sufficiently 
low magnitude do not stabilize the rotor, in fact, the rotor 
RPM that results in ground resonance is increased slightly 
over the rotor RPM where the linear rotor would experience 
ground resonance. For this reason, the Duffing replacement 
of auxilary lag damping produces less hub motion for rotors 
with a lag frequency in the upper range, 0.6 or 0.7 per rev, 
vice the lower 0.3 per rev range, which produces 


unacceptable hub and blade motion. 
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2% Comparison with Theory - Baseline Rotor 


A Coleman stability plot for the baseline rotor 
follows. This plot was produced from Matlab? code that was 


modified from Rafanello [Ref. 9]. One can see the region of 
instability between 225 RPM and 266 RPM. This plot is for 


the baseline rotor with no Duffing stiffness employed. 


lity Plot - Baseline Rotor 
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Figure 23 - Coleman plot for baseline rotor 


One can see the duplication of this region of 
instability in the simulation results, Figure 24. This 


stability plot matches the plot from Coleman quite well. The 
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region of instability begins at 225 RPM in similar fashion 
to the theory. The NPS simulation becomes stable at 261 RPM 
where the theory does not predict stability until 266 RPM. 
This difference matches the results found by both Robinson 
and Rafanello in that Coleman's linear analysis, as applied 
by Deutsch, produces slightly conservative stability 


criteria for ground resonance [Refs. 1, 9]. 
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Hub Damping As A Function Of Omega 













| Omega, /Omega-0.6, StructDamp, ,,=0.02 


lf aaa ona ono ananasa aaa alna ananasa 
$ . e . . . 
; ¿ ; i i 

Nvl—— ——À M vna 
.02 ? 7 ; ; ž 


egion of Instability 


0 A A O A O A RI A s OV 


< 3 5 : 
—————————ÓQ———————————————— rones 
: ; : : 


Hub Response Damping Ratio 


210 220 230 240 250 260 270 


Omega (RPM) 


Figure 24 - Baseline rotor region of instability 
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3 Comparison with Theory - Duffing Rotor 

The following Coleman stability plot is for the 
baseline rotor with 8 degrees lag amplitude with a Duffing 
coefficient 20 times the linear stiffness. One might note 
that Coleman's is a linear analysis. This linear analysis is 
applied to the nonlinear Duffing rotor at a single instance 
in time with a specific lead angle, and is not valid at any 
other blade lead angle. One can immediately see that the 
nominal rotor speed of 239 RPM is now stable. The Duffing 
stiffness in the rotor continuously drives the unstable 


region above the RPM of the rotor. 
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Coleman Stability Plot - Duffing Rotor 
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Figure 25 - Coleman plot for Duffing rotor 


4. The Region of Instability - Duffing Rotor 

The overall effects of including Duffing stiffness in 
the lead/lag dynamics can be more easily discerned in Figure 
26. Hub lateral stability at a variety of rotor speeds is 
plotted for four different Duffing blade stiffnesses. The 
introduction of Duffing stiffness at 10 times the linear 
stiffness greatly increases overall hub motion stability, 
but does not stabilize the rotor. A Duffing stiffness just 
over 30 times the linear stiffness would provide a stable 


baseline rotor for all rotor speeds. 
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Hub Damping As A Function Of Omega 
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Figure 26 - Change in region of instability with Duffing stiffness included 


It should be noted that in Figure 26 that the damping 
values plotted on the Y-axis are linear approximations for 
damping obtained from a nonlinear system time history (for 
all but the K,/K=0.0 line). The Hilbert transform method was 
used to determine damping from the hub lateral time history 
for each point on the lines depicted. The last 8 seconds of 
a 9 second simulation were used as inputs to the damping 


analysis. Example time histories for the Kd/K=20.0 rotor 
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follow in Figure 27. These plots were taken at the three 


Separate rotor speeds labeled. 
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Figure 27 - Example Hub Time Histories 
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5. Duffing Rotor Effect on Structural Moments 
Figure 28 depicts the blade restoring moment due to the 
combination of linear and nonlinear stiffness as a function 


of blade lag displacement. For the baseline rotor operated 


at the center of instability, the lag amplitude was 


approximately 8 degrees. This condition increases the 


restoring moment less than 50% from the linear spring only 


Case. 
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Figure 28 - Lead/lag restoring moments 
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Keep in mind, however, that a blade experiencing these 
increased loads due to root end stiffness encounters no lag 
moments from the now removed auxiliary lag damper. Also, the 
8-degree lead angles seen in the baseline case with Duffing 
stiffness (Figure 16) was for a "worst case" at the center 
of the region of instability. The center of instability is 


an unlikely nominal design RPM for a rotor, even one with 
nonlinear stiffness. For the cases shown with increased Q, 


the blade lead angles were on the order of a degree, which 
shows virtually no additional restoring moment on this plot. 
The addition of Duffing stiffness at reasonable rotor 
speeds, therefore, is not producing blade bending moments 


that current rotor designs do not already handle. 


6. Stiff-in-Plane? - Duffing Rotor 

In ground resonance analysis, one additional question 
arises: Since a stiffening term is being employed, is the 
rotor becoming stiff-in-plane? A stiff-in-plane rotor is one 
that has a lead/lag natural frequency greater than the 
rotating frequency of the rotor head. For the baseline rotor 
with a rotary spring Duffing stiffness coefficient 20 times 
the linear stiffness coefficient, the rotor only goes stiff- 
in-plane when the blade lag amplitude exceeds 6 degrees at 
239 RPM. For the baseline rotor with Duffing stiffness 


employed, the rotor can go stiff-in-plane. If the rotor is 


Su 


operated off the center of instability, however, it does not 


go stiff-in-plane. 
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Figure 29 - Soft in-plane frequencies at 239 RPM 


As depicted in Figure 30, in this case, the rotor goes 
stiff-in-plane when the blade lag amplitude exceeds 12 
degrees for this baseline rotor: with Duffing stiffness 
employed at 310 RPM. 

At the increased rotor operating speed of 310 RPM, the 
limit cycle lag amplitude was less than a degree. For lead 


angles of this small amplitude, the rotor sees very little 
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migration towards  stiff-in-plane from its initial lag 
frequency. Since the exact same baseline rotor is modeled, 
note that the lag natural frequency is changed due to the 


increase in rotor RPM. What was a 0.6 per rev lag frequency 


rotor is now a 0.6*(239/310)0 - 0.46 per rev lag frequency 
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Figure 30 - Soft in-plane frequencies at 310 RPM 
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G. CHAOS ANALYSIS 

1i. Response to Similar Initial Conditions 

Duffing's equation ls well documented In the 
literature. In Ueda [Ref. 19], parameter combinations that 
produce chaotic response from  Duffing's equation are 
summarized. Since chaotic response in the blade lag motion 
is undesired, a check of the simulation time histories must 
be completed to rule out chaotic response for the Duffing 
rotor presented in this paper. 

One characteristic of chaotic systems is vastly 
differing time histories resulting from nearly identical 
initial conditions. Figure 31 is a Poincare section of the 
lateral response for the baseline hub with  Duffing 
stiffness. Three initial conditions of hub displacement are 
plotted. The hub response depicted in this plot is not 
chaotic; the sampled points on the phase plane are nearly 


identical for the three cases. The sample frequency for this 


Poincare section is the regressing lag frequency, (Q0-0,). 
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Figure 31 - Poincare section of hub lateral motion 


2. Poincare Phase Plane and Phase Space Analyses 

In addition,  cbaotic "motzon "would “not follow a 
prescribed form in phase space. Figure 32 is the blade 
lead/lag motion plotted in 3-D phase space with the blade 
lead angle rate on the vertical axis, the lead angle on the 
radial axis, and the phase plane rotated at the lead/lag 
natural frequency. The initial conditions of zero lead angle 


and rate quickly begin to track a predictable path in the 


phase space, showing no tendency towards chaos. 
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3-D Phase Space 





Figure 32 - Phase space of blade lead motion 


H. DROOPING ROTOR STABILITY 


1% Nonlinear stiffness provides rotor stability for a 
wide range of physical parameters 


Rotor design engineers must insure rotor stability for 
a wide range of operating parameters. The fuselage rigid 


body natural frequency in roll is a strong function of the 
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fuselage inertial properties. Different fuel loads, cargo 
and passenger configurations effect fuselage rigid body 
natural frequencies. Landing gear properties, such as strut 
cleanliness, strut oleo servicing, and tire inflation can 
also effect these frequencies. The engineer must work to 
obtain a stable rotor under all these conditions. In 
addition to fuselage characteristics changing, rotor blade 
properties can change. Helicopter operators can change trim 
tab and pitch link settings, effecting the way the blade 
tracks. They can also change blade tip weights in the 
dynamic balancing process, which has a strong effect on 
blade inertias. 

Historically, larger lead/lag dampers than were 
necessary were employed in final designs in an attempt to 
Maintain rotor stability over all possible operating 
conditions. The dampers were used as an engineering "safety 
net." In Figure 33, the Duffing stiffness in the rotor is 
used as a similar safety net for rotor speed reduction. 
Reduced rotor speeds could be experienced as a result of 
system failure or during maintenance checks, in addition to 
rotor engagement and disengagement. This figure compares a 
motor With Duffing stiffness BO a rotor with linear 
stiffness only. As rotor RPM is reduced, both rotors enter 
the region of instability and experience resonance. The 


linear rotor quickly reaches hub displacement amplitudes 
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that could cause catastrophic failure. The Duffing rotor 
encounters smaller hub displacements, avoiding loss of the 
helicopter. In this case, the Duffing stiffness is used as a 
rotor safety feature, avoiding the large hub displacements 


that result in failure as the rotor speed is reduced. 
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Figure 33 - Reduced rotor speed comparison 
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IV. LINEARLY LINKED ROTOR SYSTEMS 


A. INTRODUCTION 

The helicopter in ground resonance could be compared to 
a simple harmonic oscillator. The rotor acts as the exciting 
force, and the rigid body motion of the helicopter is the 
response. If one desires to avoid resonance, one could take 
one of two approaches to the problem. Frequencies of the 
system could be modified, which would avoid coalescence. 
Either the exciting force frequency (the rotor regressing 
mode) or the oscillator natural frequency (the fuselage 
rigid body natural frequencies) could be changed to avoid 
driving the system into resonance. 

The other approach that could be used is the 
modification of the resonant mode shape. Changing the mode 
shape naturally changes the mode’s resonant frequency. By 
linking the blades together, the mode shapes of the blades’ 
lead/lag motion are changed. Modifying the mode shape, in 
turn, alters the blade natural frequency in lead/lag and the 
regressing mode of the rotor. Therefore, by linking the 
blades, the potential exists to avoid ground/air resonance 
by detuning the rotor-body dynamics, thus avoiding ground 
and air resonance in a similar manner to the nonlinear 


stiffness approach already discussed. 
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B. ROTOR WITH LINKED BLADES 

A four bladed rotor simulation was developed similarly 
to the  Coleman-like three bladed rotor used in the 
damperless rotor investigations and the Comanche five bladed 
rotor used in the second chapter. This baseline rotor was 
intentionally chosen to be unstable. The eguations of motion 


for this new four bladed simulation were also derived using 
Lagrange's method in a Maple? worksheet. Unlike the other 


rotor simulations, this new four bladed rotor had the 
additional capability of linking the blades, either 
elastically or rigidity. In addition to linear stiffness acma 
function of the angle made between the blade and its nominal 
position with respect to the hub, in this case, rotary 
stiffness terms were added as functions of the lead/lag 
angles made between adjacent blades. All linked stiffness 
and damping values were linear. 


If traditional linear root-end stiffness were described 
as M=kr*C, for the n" blade, then this new interlink 
stiffness could be described as: 

Mo. ce) 
where Kae Ís the interlink stiffness (moment/angle), and 


(G., -G.) and  (6,-G) are the instantaneous lag angle 


differences between the nth blade and its two adjacent 


neighboring blades, (n+1) and (n-1). The simulation was also 
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given the capability of linking the blades with viscous 
dampers which were functions of blade lead angle rate 


differences: 


M, = Cinern+1) CR x s NEC nera) Gis -¢,) 

Investigations in the potential payoffs of inter-blade 
linking follow in the next three sections. The baseline case 
properties for the four bladed interlink rotor may be found 


in the following table: 


Table 3 
Q = 25.0 rad/s M 5 160.0 s1 M aaae = 4.0 s1 
(239 RPM) (5150 lbs) (129 lbs) 


Linear Damping in Lag = 1.0% 


Linear damping in Hub Lateral Motion = 2% 





Time histories for the baseline case blade lead angles 
may be found in Figure 34. This case also has been purposely 


set inside the region of instability, and clearly diverges. 
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Figure 34 - Four bladed interlink baseline 


C. LINKED BLADE PARAMETRIC STUDIES: 


1% Linked blade & hub response varying blade 
interlink stiffness for baseline rotor 


The effects of varying blade interlink stiffness for 
four different hub damping values can be seen in Figure 35. 
For the cases presented, each blade is elastically linked to 
each of the two neighboring blades in series. For very low 
malen Nini MetifinesS v RO (K, as low, the simulations ane 
unstable, similar to the baseline case, as would be 


expected. As the interlink stiffness is increased, however, 
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a sharp increase in rotor stability occurs in the region 
where the blade interlink stiffness is of the same order of 
magnitude as the flexbeam chord-wise stiffness, K,,,,.. For 
these cases, cereale unstable rotors were stabilized for 


all but the lowest hub damping of 1%. 


Blade Response Damping vs Blade Link Stiffness Ratio 
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Figure 35 - Blade elastic interlink results 
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2 Linked blade & hub response varying blade 
interlink damping for baseline rotor 


Interlinking  dampers was also attempted. All four 
blades were interconnected by similar dampers, thus putting 
the interlink dampers in series. Recall for the baseline 
four bladed rotor that linear blade damping was only 1.0$. 
For inter-blade damping, the same four hub values are 
plotted in Figure 36. The interlinking of blade dampers also 
proves effective. Unlike the case of the  interlinked 
stiffness, this increase in stability is to be expected 
since the overall system damping has been increased, even 
though the additional dampers are implemented in less than 


traditional means. 
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Figure 36 - Blade damper interlink results 


3. Linked blade & hub response with uneven interlink 
stiffness for baseline rotor 


Looking at  unsymmetrical blade interlinking was 
accomplished so as to change the relative phasing of the 
blades in ground resonance. Note that for the divergent 
baseline case (no linking), each blade moves 90 degrees out 
of phase with its neighbor in the case of a four-bladed 


rotor. In general for an n-bladed rotor, Coleman showed that 
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during divergence the blades moved at a phase angle of 2n/n 


(radians) or 360/n (degrees) with respect to each other. It 
was expected that disrupting this phasing would change the 
magnitude of the rotor CG offset from the center of 
rotation, in turn reducing the coupling into the fuselage 
degrees of freedom. This is due to the fact that we can 
change the natural frequency of a system by altering its 
eigenvalue ee, the spring rate) or by altering its 
eigenvalue (changing the mode shape). In Figure 37, the #1 
and #2 blades alone are linked. This action does reduce the 
amplitude of the blade lead angles over the baseline 
unstable case, but all cases using this method diverged, 


nonetheless. 
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Figure 37 - Blades +1 and +2 linked 


Linking the blades in pairs, however, showed far 
different results. In Figure 38, the #1 and #4 blades are 
linked, and the #2 and #3 blades are linked. This 
interlinking in pairs stabilizes the rotor. For a linearly 
unstable rotor, blade limit cycle amplitude is only 2 


degrees. 
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Blades Four/One and Two/Three Linked 
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Figure 38 - Blades #4&1 and +283 linked 


By linking the blades in adjacent pairs, the four 
bladed rotor acts more like a two bladed rotor, unable to 
enter ground resonance. This method of stabilizing the rotor 
shows promise in that the aerodynamic and handling response 
characteristics of a multi-blade rotor can be maintained, 
but the rotor exhibits the mechanical stability normally 


found in two bladed rotors. 
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4. Interlink Dampers Versus Conventional Dampers 

In Figure 39, traditional blade dampers are compared to 
interlink dampers. Essentially, the four lead/lag dampers 
that might be found on a rotor have been detached from the 
hub and remounted between the four blade roots. Identical 


damping values were used for the two lines plotted. 


Blade Response vs Auxiliary Damper Strength 
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Figure 39 - Conventional dampers versus interlink dampers 


The interlink damping proves more effective for a given 
damping value. Instead of the dampers acting independently 


for each blade, the new system takes advantage of the fact 
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that the blades are 90 degrees out of phase with each other, 
essentially amplifying the input to the dampers. Since 
weight is always a driving issue in aircraft design, if the 
dampers were mounted between the blades, instead of from the 
blade to the hub, smaller dampers could be installed and 
still maintain the desired damping in the blade chord-wise 


dynamics. 
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V. ROTOR WITH UNEVEN BLADE SPACING 

As long as unsymmetrical rotors are being investigated, 
the effect of modifying blade spacing was also attempted. 
Blade motion time histories from a 60-120-60-120 degree 
spread rotor can be found in Figure 40. This rotor layout is 
Similar to the AH-64 Apache tail rotor, where the four 
blades are not positioned 90 degrees apart, but instead 
arranged in an X pattern, not a perfect cross. The Apache- 
like rotor also exhibited blade divergence. The blade pairs 
diverged at different rates, but the same overall response 
was noted: mechanical instability was not avoided for any of 
the rotor speeds investigated for the case where uneven 


blade spacing was introduced. 
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VI. GROUND AND AIR RESONANCE ANIMATION 
Obviously, numerous cases and configurations were 


simulated in the course of this research. All these 


simulations were performed on a PC running Simulink?. 


Utilizing the numerical analysis capabilities of Matlab? 


including the Hilbert damping analysis, numerous 
quantitative parameters were generated. 

In an attempt to provide a better qualitative picture 
for what the rotor was doing, an animation routine was 
written that would show the engineer visually what the 
blades and hub were doing in the horizontal plane. The 
animation routine written as part of this research plots the 
rotor blades as simple lines at the computed lead angle from 
the simulation at every time step. The depicted rotor does 
not rotate, so that blade action can be observed. The center 
of rotation is also changed. The animated blade center is 
plotted at the hub X and Y position with respect to the 
inertial frame, also obtained from the simulation at every 
time step. In this way, the outward spiral of the hub is 
also depicted as part of the animation. 

Though no numerical results are obtained, the animation 
routine was found very useful in providing an overview of 
rotor states. Four example frames of this simulation are 


shown in Figure 41. These frames show the progression of 
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blade movement in a typical ground resonance scenario. 
Notice how the animation illustrates the spiral motion of 
the CG of the rotor around the center of rotation. The 90 


degree blade phasing is also easily discernible. 
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Rotor Simulation Animation 


Simulation time: 
9.9167 sec 


à * t— 60? in Blade "1 Lead Cycle 
(blade "1 is at 6 o'clock in figure) 


9.9833 sec 


o *t- 120? in Blade *1 Lead Cycle 


10.1167 sec 


œ * t = 240? in Blade "1 Lead Cycle 


10.1833 sec 
@ * t= 300° in Blade “1 Lead Cycle 





Figure 41 - Rotor animation 
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VII. DAMPERLESS ROTOR DEVELOPMENT CONCLUSIONS 

With the advent of hingeless and bearingless main rotor 
helicopters, interest in nonlinear rotor dynamics has 
increased significantly in recent years as an area to 
research for improved rotor stability. With the rapid 
evolution of high-speed processors and recent developments 
in software technology, we now have the tools available to 
explore rotor system nonlinear dynamics and potentially 
achieve helicopter rotor mechanical stability. 

This dissertation furthered the development of a non- 
linear analysis tool that has the capability of modeling a 
rotor's nonlinear material properties. Most importantly, the 
simulation did not incorporate assumptions or ordering 
schemes that limit the number of terms in the equations of 
motion and subsequently skew the output. In the simulation 
used, not a single term was eliminated. 

Potential methods for passively controlling lead/lag 
motion of helicopter blades have been developed. One of the 
methods presented utilizes Duffing-type nonlinear stiffness 
in the root end of the rotor blade in addition to linear 
spring and damping properties. The Duffing spring provides a 
restoring force that is a cubic function of the blade 
lead/lag deflection. This nonlinearity shifts the first in- 
plane natural frequency of the blade as a function of the 


blade lag angle. This shift in blade natural frequency 
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effectively varies the regressing mode frequency. The 
resulting shift in the regressing mode prevents coalescence 
with hub resonant frequencies. For a linearly unstable 
rotor-fuselage system, inclusion of the Duffing spring 
results in limit cycle behavior in the blade lead/lag 
degrees of freedom and in hub translational motion. 

An analysis of the limit cycle motion effect on hub 
vibrations was also completed. The replacement of blade 
auxiliary lag dampers with nonlinear springs results in 
increased hub vibration levels if the resulting system is 
highly unstable in the linear analysis. Rotors that are 
marginally stable, or stable, in the linear analysis produce 
vibration levels that are not uncommon in production 
rororeraft. 

The primary conclusions of this research on damperless 
E@EOrS dre: 

e Application of nonlinear stiffness properties to the 
blade lead/lag degrees of freedom can effectively improve 
helicopter ground/air resonance stability 
characteristics. 

e Stability characteristics are improved through frequency 
mismatch, not by dominating system dynamics with high 
damping. 

e Blades with low lead/lag natural frequencies have greater 


potential for producing unacceptable vibration levels in 
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both dynamic and non-rotating components if Duffing 
stiffness is employed. 

Interblade rigid and elastic coupling has the potential 
to eliminate helicopter unbounded blade motion. 

In this case, the interblade coupling achieves stability 
by altering the phase of relative blade response. That 
is, the response eigenvector has been altered. 

Interblade viscous damping, selectively applied, can be a 
more efficient use of blade auxiliary dampers than the 


conventional blade-to-hub arrangement. 


Uneven blade spacing, for the case examined (60-120-60- 
120), did not prove effective in avoiding helicopter 


ground and air resonance. 
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APPENDIX A. COMANCHE FROUDE SCALE MODEL SIMULATION 


Comanche Rotor Design Overview 


Engineering data on the Boeing/Sikorsky RAH-66 Comanche 
Froude scale wind tunnel model was obtained from Sikorsky 
Aircraft as part of a research program initiated through the 
National Rotorcraft Technology Center (NRTC) [Refs. A1, A2]. 
Additional material on the Comanche rotor is provided by 
Tarzanin and Panda [Ref. A3]. 

The Boeing/Sikorsky report obtained by NPS outlines 
testing procedures, results, and basic rotor design geometry 
of the Comanche model rotor. The Comanche rotor design 
utilizes a bearingless main rotor (BMR) flexbeam arrangement 
for blade lead/lag, flap, and feathering degrees of freedom. 
A cuff, which is stiff in torsion and bending, encases the 
flexbeam and provides input torques in the feathering axis 
to the blade at the outboard end of the flexbeam, at the 
cuff/blade joint. A damping device called the snubber-damper 
is mounted between the inboard end of the cuff and the hub 
to provide additional damping and stiffness in the blade 
lead/lag degrees of freedom. A depiction of this rotor is 


given in Figure Al, reproduced from Panda. 
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Secton 
Ffuidiastic? Configuration 





Figure A1 - Froude scale Comanche rotor design [From Ref. A3] 


Description of Nonlinear Snubber-Damper 

The snubber-damper is made up of alternating layers of 
metal and elastomer, which provides damping when the layers 
are put in shear. Unfortunately, the damping and stiffness 
values for the snubber-damper are not constant; they are a 
strong function of the snubber-damper displacement 
amplitude. The initial snubber-damper designs provided very 


little damping for small deflections. 
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In addition to the low damping values, snubber-damper 
stiffness tended to rise sharply for very small amplitudes. 
Plots of the snubber-damper stiffness and loss factor were 
published by Panda and follow as Figure A2, [Ref. A3]. The 
"Elastomeric Damper" depicted is the device used in the 


1992-1993 wind tunnel tests and the "Fluidlastic?" damper 


was the result of a redesign, tested in 1995. 
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Elastomeric Damper 


Fluidlastic® Damper 


K Stifíness (Ib/In) 


0.01 0.02 0.03 0.04 0.05 0.06 


Displacement Amplitude at 10 Hz (in) 


7 Fluidlasti? Damper 


Qo mm CD CUM Gem mono mev sms 


Loss Factor 


Elastomeric Damper 


0.01 0.02 0.03 0.04 





Figure A2 - Snubber-damper stiffness and loss factor [From Ref. A3] 


Inherent in the design of these devices, the snubber- 
dampers exhibit hysteresis in their blade restoring force as 


a function of displacement. Panda reported the Comanche 
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elastomeric  snubber-damper as having sufficiently low 
damping at small displacements that this device caused limit 
cycle oscillations in the blade chord-wise motion. This 
limit cycle motion produced undesired vibrations in the 
coupled rotor-fuselage system. The Fluidlastic® snubber- 
damper alleviated the severity of the limit cycle motion 
problem. One can see less reduction in the Fluidlastic® 


damper loss factor at small displacement amplitudes in the 


lower plot of the previous Figure. 


Flexbeam and Cuff Modeling 

In order to accurately model the 1/6 scale Comanche 
wind tunnel test rotor, flexibility in the blade and 
flexbeam had to be accurately modeled. Initially, a linear 
beam approximation was used. This was done primarily as a 
first iteration in the model design process. The linear 


theory also provided data for verification of a Myklestad 
method that was coded in Matlab? as an additional task in 


this research. 
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Comanche Model Assembly 
The Comanche 1/6 scale flexbeam and blade were modeled 
using the Myklestad program. Example output of a l-inch 


blade tip deflection in lag follows as Figure A3. 
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Figure A3 - Combined flexbeam/blade analysis 


Using the Myklestad program, the equivalent offset was 
found to be a constant value of 8.52 inches from the center 
of rotation, with the snubber-damper not installed. This 


equivalent offset is highlighted in Figure A4 by a small 


circle. 
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Figure A4 - Flexbeam/blade deflection 


From Figure A4, it can be easily seen that the outboard 
aerodynamic portion of the blade moves as a rigid body for 
in-plane bending at the frequencies of interest. Therefore, 
the portion of the blade outboard of the equivalent offset 
was modeled as a single degree of freedom. Since the slope 
change in the combined model is at the root end in the soft 
region of the flexbeam, an articulated type attachment could 


be used for mathematical modeling of the hingeless rotor in 


the Simulink® simulation. 
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Using the blade mass distribution data provided, a 
blade mass and moment arm were computed. A blade moment of 
inertia was calculated from this data about the equivalent 
offset. This inertia was then broken down into an equivalent 
mass, m,, and moment arm, 1l. The blade was modeled as a 
point mass at a distance 1 from the offset, such that the 
total mass of the actual blade matched the point mass value 
in the model. For accurate modeling of the  lead/lag 
dynamics, it was desired that e*S/I for the actual blade 
match e/l for the point mass model, where S is the first 
mass moment of the blade and I eguals the second mass moment 
of the blade. The blade model was determined to be a 0.0129 
slug point mass at a radial location of 30.06 inches. 

Model parameters for the portion of the blade outboard 
of the effective offset are compared to the actual blade in 


the following table: 


Table Al 










NPS Model 
Blade/cuff moment of inertia, I 
Blade/cuff mass moment, S 
Blade/cuff mass, m, 


TAO 2135407 1n 


Blade lag offset "OHA sl 
im oc ni 











8.5435 XEM 
30.0642 in 






Point mass radial location 





The cuff is far stiffer than the flexbeam, providing an 
essentially rigid sleeve around the flexbeam. For lead/lag 
motion, the flexbeam simply bends inside the interior 


dimensions of the cuff. Since the snubber-damper attaches at 
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the inboard end of the cuff, the snubber-damper 
displacements were considered to be a function of an 
extension of the cuff/blade joint angular displacement and 


off-axis deflection, as seen below in Figure A5: 


Joint angular displacement 


Cuff-blade joint e 


Joint off-axis deflection 


Snubber deflection 





Figure AS - Snubber-damper motion as a function of cuff/blade joint motion 


Assuming a constant equivalent offset for blade 
lead/lag motion, a geometric gain was computed for snubber- 
damper displacement as a function of blade lead/lag angle. 
Using the physical properties of the blade root end, the 
geometric gain was computed to be 5.937 inches of snubber- 
damper deflection in shear for every radian of blade lead, 


assuming small lead/lag angles. 


Snubber-Damper Model 


ERR 


The NPS nonlinear rotor simulation was modified to 
incorporate nonlinear snubber-damper coefficients in both 
stiffness and damping. The loss factors made available in 
the Boeing/Sikorsky report were first converted to 
dimensional damping coefficients. The loss factor is defined 
as the out-of-phase force divided by the in-phase force. 
Since the experimental loss factor data in the report was 
recorded at a fixed driving frequency of 10 Hz, the loss 
factor data could be converted to dimensional damping 


values. The dimensional damping value may be expressed as: 


pe 
w 


where 7 is the loss factor, k is the nonlinear snubber- 


damper stiffness, ® is the forcing function frequency, and 


x is the snubber-damper displacement. 

Traditional hydraulic and viscous damping methods were 
attempted in the nonlinear  snubber-damper mathematical 
model, ain addition to higher order curve fits. A fourth 
order fit, Cc(x)=c,X'+C,X'+C,X +C,x+C,, in displacement was 
decided upon as the best trade-off between accuracy and 
model order reduction. A linear relationship was used for 
snubber-damper displacements greater than 0.6 inches. 

The polynomials formulated for both the snubber-damper 


damping and stiffness formulations follow: 
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Table A2 


1992 Elastomeric snubber-damper: 
c=-842875x +87576. 8x -2847 .07x°+25.01x+0.47 
k=336241084x -38666827x +1547897x -26483.84x-233.32 


1995 Fluidlastic® snubber-damper: 
cz116695.68x -14022.78x 1609.0x-12.05x10.59 
k=38512436.84x -4491019.58x"+183336.70x"-3147.25x+79.71 


Where ec is in lbs/zn/s, k in lbs/in, and xen inches. The 














stiffness and damping results are compared to the snubber- 


damper laboratory data provided in the plots that follow: 
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Figure A6 - Snubber-damper stiffness and damping properties 
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NPS Simulation Results for the 1/6 Scale Comanche 
Example results of the NPS rotor simulation follow. The 


inputs for these plots may be found in the following table: 


Table A3 - Helicopter physical and aerodynamic parameters 


Omega=870* (pi/30) 
el=8.5235 

1=21.5407 

Rs=2.5984 

No lag stops modeled 

















mb-0.0129 
$ Effective mass of fuselage, slugs | Mx- 0.6452 

My= 1.3802 
Phi(1:5059[0:2:74' 6 8]*BHE 


$ geometric gain of snubber gg=5.1121 
movement from lag, in/rad 

$ Spring stiffness polynomial for 
lead-lag, in*lbs/radian) (1995 
snubber) 













| 


eam=[0,0,0,0,660.9141] 
Ksnub=[38512436.84*gg”4, 
-4491019.58*99°3, 183336. 79*9905 2, 
-3147.25°59,.79. 717 

*(e1-Rs)^2 

Kpoly=Kbeam+Ksnub 


% Lead-lag stop spring Ks=0 
constants, in*lbs/radian No lag stops modeled 


% Damping polynomial in lead-lag Cpoly=[116695.68*gg"4, 
(in*lb/(rad/sec)) (1995 snubber) -14022.78*g9"3,609.0*gg"2, 
-12.05"gg,0.59] "(el-Rs)"2 


Linear springs in translation, K(1)2297.38 
WOSAN K(2)2217.93 


Linear damping in translation, EE 
(1bs/(in/sec)) c(2)=1.88 


Non-linear damping in translation, v=0 
lbs/(in/sec)^2) 


Fuselage initial displacements, in E EH 
(all other IC's = 0.0) xY1-0.1 









































The fourth order snubber-damper model is employed in 
the simulation depicted in the following time histories. All 
blades are given initial conditions of zero lead angle 
displacement and zero velocity. All rotor and hub dynamics 


are initiated by a 0.1 foot displacement in both the x and y 
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dimensions of hub motion. Hub response from these initial 
CONG TION m houwn In Pigure A7. This initial Conditien is 
equivalent to a lateral restraining force on the hub of 
36.87 lbs. being released at time t-0.0 in the simulations. 


This comparatively large impulse was set to insure coupling 


of the rotor/body dynamics. 


Hub Translational Time Histories 


| | | — Lateral 
0.08 WEE Longitudinal d 


0.1 





Displacement (in) 


Time (s) 


Figure A7 - Hub translational time histories 


The Comanche wind tunnel tests utilized a shaker in the 
fuselage tO initiate rotor/body coupling prior to recording 


time histories for analysis. Since the initial conditions 
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varied for each test in the tunnel, duplication of the exact 


initial conditions was impossible with the simulation, and 


the fixed initial conditions on the hub displacement were 


used. The Boeing/Sikorsky report describes this process as 


EO Ows: 


"The blade lead-lag mode was excited by the 
swashplate cyclic shaker operating at the fixed 
system lead-lag regressing frequency. When 
required, a fine adjustment was made to maximize 
the lead-lag response of the model. The amplitude 
was monitored (on line) by a spectrum analyzer. 
When the forced response was judged to be maximum, 
the shaking was stopped, marking the beginning of 
the transient decay. The following 10 seconds of 
chord bending and fuselage roll-pitch data was 
collected and analyzed by a moving-block method to 
determine the frequency and damping ratio of the 


mode being investigated." [Ref. A2] 


Response to the relativey stiff rotor blades is shown 


in Figure A8. The Comanche lead/lag frequency was reported 


to be 0.710, for the 1995 system, 0.6920, for the 1992 


System. 


S] 





0:1 
: i : : : : Blade 1 
0. 08 D NL A KISS sta en a F Aa ds Ue esie plici: Blade 2 n 
Blade 3 
0.06 tika 2 Blade 4 | 
r § | ll Blade 5 
O) 0.04 FIBI - 
Z 
~ 0.02} 
© E 
O 4 
c “ 
za -0.02 
Ke . D 
D 
J -0.04 
-0.06 
-0.08 E eeben 
-0.1 
0 1 2 3 4 5 6 7 8 9 
Time (s) 


Figure AS - Blade lead angle time histories 


Damping for each of the five blades was determined 
using a Hilbert transform method [Ref. A5]. The damping 
value for the blades was then averaged for the final 8 
seconds of the 9 second time history. This average was then 
converted to a fixed system damping value replicating the 
method used in the Comanche test: 

"The rotating system damping ratio is converted to 

a fixed-system damping ratio by multiplying by the 


ratio of the rotating-system lead-lag mode frequency 
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(0) to the fixed-system lead-lag regressing mode 


frequency (Q-0x)." [Ref. A2] 


An example output of the Hilbert damping analysis for 


the #5 blade follows: 


Hilbert Method Damping Determination 
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Figure A9 - Blade motion damping analysis 


The top subplot is the time history of blade lead angle 
that was analyzed. The second subplot is the natural 


logarthm of the absloute value of the Hilbert transform of 


do 


the signal. A first order fit to this plot is also shown. 
The damping of the signal was then computed by dividing the 
negative of the slope by the time history frequency. The 
primary frequency of the time history is determined by power 


spectral density using Matlab?'s discrete fourier transform 


function, "fft.m." The results of the power spectral density 
are shown in the third subplot. This Hilbert transform 
method works well for signals dominated by a single mode. 
Smith obtained more accurate results with the Hilbert 
transform method than with Moving Block for signals with a 
single mode present [Ref. A5]. Results of the Hilbert 
transform damping analysis of the two hub degrees of freedom 


follow in Figures A10 and A11. 
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Figure A10 - Hub X direction damping 
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Figure A11 - Hub Y direction damping 


The frequencies obtained from the NPS simulation are 


summarized in the following table: 
Table A4 
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Comparison with Test Data 

The goal of the application of the NPS simulation to 
the Comanche scale rotor data was the matching of fixed 
system damping values to the test data and to UMARC. The 
plot that follows reproduces the Boeing/Sikorsky report's 
Figure 50b, from the 1992 test results [Ref. A1]. The 1992 
test data was chosen as a more rigorous test of the NPS 
simulation, since the snubber-damper used in the 1992 tests 
exhibited greater nonlinear behavior. The first plot, Figure 


A12, shows the NPS simulation's initial results: 
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Replication of Boeing Sikorsky Figure 50b, 1992 
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Figure A12 - Initial NPS simulation results 
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The geometric gain, the snubber-damper shear to blade 
lead angle ratio, was considered a soft parameter and varied 
in an atttempt to obtain more accurate matching with the 
existing data. 

As the geometric gain was increased, slightly better 
results were obtained. The following plot shows the results 
of increasing the geometric gain by a factor of 1.5. From 
these results, the wind tunnel model, it is assumed, was 
Seeing more snubber-damper movement than the modeling method 
assumptions had called for. This soft parameter, the 
geometric gain, is related to the multiple load path problem 
alluded to earlier. In brief, the basic assumptions 
underestimated the amount of snubber-damper motion with lead 
angle changes. Recall only lead/lag motion is modeled, and 
the bearingless Comanch rotor model experienced coupling in 
the lead/lag dynamics from other model degrees of freedom. 

With the advantage of having the test data results a 
priori, the model was tuned to best replicate the wind 
tunnel data using the geometric gain as a tunable parameter. 
With the higher geometric gain mentioned, the NPS simulation 
appeared to replicate the model e slightly better than 
UMARC at the lower RPM's tested. At higher rotor speeds, the 


difference between the two simulations is negligible. 
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Figure A13 - Final NPS simulation results with increased geometric gain 


Comanche Modeling Conclusions 
e If the physical rotor exhibits migration of the 
effective offset with varying lead angles and 
velocities, forcing a hingeless rotor into an 
articulated model contaminates the results. 
e The inclusion of aerodynamics, coupling in the rotor 
lead/lag, flap, and torsion should be included in 


the modeling of a bearingless rotor. 
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e Modeling the blade portion that is outboard of the 
effective offset with only two parameters, blade 
length and tip mass, is incapable of matching all 
the inertial properties required to accurately 


simulate both lead/lag and hub motion. 


e Additional degrees of freedom are required in the 
lead/lag dynamics to better approximate multiple 


load paths from the blade to the hub. 


e For a fledgling code, the NPS simulation shows 
potential in modeling rotors with nonlinear physical 
parameters and further development should be 


considered. 
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APPENDIX B. MYKLESTAD BEAM MODELING METHOD 


In 1944, Myklestad published a method for finite 
element analysis of beams in bending for vibrations 
applications in his book, "Vibration Analysis" [Ref B1]. His 
method breaks a beam into segments that have essentially 
linear behavior over the segment length. Equations of motion 
of each of the separate elements are then tied together as 
separate degrees of freedom. In Myklestad's method, each 
beam element is defined as a point mass and a massless 
elastic element in bending, as depicted in Figure Bl, 


modified from Gerstenberger and Wood [Ref. B2]. 





Figure B1 - Typical Myklestad beam element [After Ref. B2] 


By defining element members in this way, a beam with 
inertial and mass properties that change as a function of 
beam longitudinal position may be well approximated. Shear, 
moment, angular displacement and beam deflection may be 


described for each element as follows: 


"A SS, +m,y,@° 
Ma: =M, + Pd ye 
Up n+l ya T 
Or, -0, Te M = 
E, JER 
2 l B 


n,n+l 


n+1 


=y + O F- M n- 
Y n+1 Yn nna yn 2E,I, n+l 3E I, n+l 


where @ is the first beam natural frequency in bending, 1,,, 
is the segment length, E, is the element stiffness, and I, is 
the element second area moment. 

Myklestad's method may be modified from the simple beam 
in bending to rotating beams by the inclusion of centrifugal 
and blade tension terms. For this case, the x-y plane is 
transformed to a radial and off-axis displacement frame, as 
dene. in Figure B2, also reproduced from Gerstenberger 


and Wood. 
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Figure B2 - Typical rotating Myklestad beam element [From Ref. B2] 


The equations for a Myklestad beam element in the 


rotating frame follow: 


C 


=S 


n n+1 


P =T Frm âT, 
+S ] 


M, = kal n+1"n, n+l 7T Ga I x,) 


ie OF d Ja nn+l ) T MV nl u 


+m, (0? +0? )x, 


Xa Xl no One n+l + TS n+l M, ea = DS ed 
where V, n+l ET. es n+l Fl nn VE and Cm n+l sl. ee S 


A program written in the Matlab? application language 


that analyzes a beam using Myklestad's method may be found 


elsewhere in the Appendicies. This program analyzes beam or 


TET 


rotor blade inertial, geometric, and stiffness data for any 


number of beam segments, up to the limitations of the 
Matlab? host program used. Both rotating and non-rotating 


beams can be analyzed. 

The program also has the capability of applying 
different boundary conditions at the ends of the beam. An 
articulated blade requires zero moment and displacement at 
the lag hinge where a hingeless rotor requires zero angular 
displacement and deflection at the root end. In both cases, 
shear and moment must be zero at the blade tip. 

Linear beam theory was utilized as a means of verifying 
output of the Myklestad bending program. In the case of the 


flexbeam, where clamped end restoring moments were needed, 


M EI 
the rotary spring constant of my was used to check the 


o 

Myklestad output. This rotary spring constant for the blade, 
namely the flexbeam root-end moment per unit tip angular 
deflection, was a required input for the NPS roter 
simulation. 

Frequencies of oscillation were checked against Den 
Hartog's text "Mechanical Vibrations" (Ref. Bl]. For 
example, the first lag EE of the hingeless blade 


determined by the Myklestad program was checked with the 


E 
cantilevered beam frequency from Den Hartog, 023.52 uc 
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The plots in Figure B3 shows typical output of the 
Myklestad program. The third subplot, the moment as a 
function of beam longitudinal position, shows a comparison 
to linear theory. The example beam used in this comparison 
has a length of 10 inches, a weight distribution of 0.1 
lb/in, and a constant EI of 10000 lb*in. It can be seen 
that the Myklestad program correctly applies the four 
boundary conditions of zero slope and deflection at the 
root, and zero shear and moment at the tip. It should be 
noted that a small correction is made to the shear at the 
tip of the beam. This is done due to the final Myklestad 
element having a net shear at the tip one beam element away 
from the final mass element. This tip correction to shear is 


common practice in applications of Myklestad's method. 
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Figure B3 - Myklestad analysis output 


The following table compares a non-rotating 
cantilevered beam to Den Hartog's closed form equations, 


verifying the frequencies obtained from the Myklestad 


program: 


Table B1 


Non-rotating cantilevered beam 69.6464 Hz 69.2986 Hz 


The Myklestad computer code follows in the appendices. 
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APPENDIX C - MYKLESTAD METHOD PROGRAM IN MATLAB 


function myklestad 


clear 
format compact 


$ Required inputs: 


Ohm=0; % Rotor rad/sec 
5olbs0.l; % tolerance in acceptable assumed frequency 
(rad/sec) 


load h66blade.mat -ascii 
el=2.5984; 

% snubber radial location, in 
radius=38.976; 

% blade radius, in 


rb-fliplr(h66blade( : , 1) " radius) ; $ in 
wb-fliplr(h66blade ( : , 2) !" ) ; $3 05782 
Elb=fliplr (h66blade(:,3)'); % lb/in 
urerb(i:17); 

wf-wb(1:17); 


BEREZEIDUI:TT):. 

figure(4) plot (r£,ELE/1000) , grid 

xlabel('Radius (in)','FontSize',16) 

ylabel('EI, in^2Klb','FontSize',16) 

title('EI vs Radial Location for Flexbeam','FontSize',16) 
muf=(wf/32.2); 

nfe-length(rf); 


Beutbd=eldiff(rE),0]; 

eanbd=e[0,dIff(rf)]; 

METfF-sum( (dinbd+doutbd)/2.*EIf)/(rf(nfe)-rf(1)); 

disp((' The flexbeam length is: ',num2str(rf(nfe)-rf(1)),' in']) 
disp([' The mean EI for the flexbeam is: ',num2str(mEIf),' lb*in^2']) 


% obtain beam natural frequency and root-end transfer matrix 
% Finds natural frequency and dynamics matrix using 
% Myklestad method for beam bending 


n=O; 


if nfe--length(EIf) 

error ('The EI vector is not the same size as the r vector.') 
elseif nfe~=length(muf) 

error(’The mass/length vector is not the same size as the r 
vector.’) 
end 


omega_1=3.52*sqrt(mean(EIf) /(mean(muf) *(rf(length(rf))-rf(1))*4)); 


% rad/sec 
omega-[0.9 1.1]*omega 1; 
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SS X ke ke ke e kv à kx kx kx kx kx X k KIHILE LOOP * * * * * * * * * x kk k kk xk kx kx * 


while (abs (omega(1)-omega(2))>tolb) 
k=k+1; 


4 Th tip = 1, all others zero 
trips 0.9 280] 
state-mykloop(tip,omega(1),rf,muf,EIf,Ohm); 


$ compute b components of dynamics matrix 
bTh=state(1,1); bX=state(2,1); bM=state(3,1); bS=state(4,1); 
$b-[bTh, bX, bM] 
% debug 
*+plotstate(2,r, state, omega (1)) pause 
% debug 


$ X tip - 1, all others zero 
e1p- (0, 7,2020), 
state=mykloop (tip, omega(1),rf,muf,EI£,Ohm) ; 


% compute a components of dynamics matrix 
aTh=state(1,1); aX=state(2,1); aM=state(3,1); aS=state(4,1); 
$a-[aTh, aX, aM] 


$ debug 
%plotstate(3,r,state,omega(1)) ‚pause 
% debug 
% compute determinant 
dynam- [bTh,aTh;bX,aX]; $ cantilevered blade 
oo EL 


omega=fliplr (omega) ; 
err=det (dynam) ; 
else 
err-[det (dynam),err]; 
% use linear interpolation for next guess at omega 
omega=([max (0, omega (1) -err(1) * ((omega(2)-omega(1))/(err(2)- 


err(1)))), omega]; 
% [spline(err,omega,0), omega]; % optional spline 
method 
end 


end % while loop 


Zeek kkk KKK KARA KKK KWH TLE LOOP * * * *** ** a A kon kx kx kx X 


% debugging 

disp([” ',num2str(k),' iterations in obtaining root transfer matrix']) 
figure(3) plot (omega (2:k+1), err,omega(2:k+1), err,'o') 
xlabel('omega'),ylabel('determinant'),grid,title('Convergence Plot') 


Mj=-bM; 
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KjzaM/rf (nfe); 


disp(' ') 

disp(’ For Simulation Inputs, the Flexbeam Produces’ ) 

disp([’ the following values for Omega = *,num2str(Ohm*30/pi),' 
RPM' ]) 

disp([’ Mj = *,mum2str(Mj),’ in*lb/rad deflection at the jornt" | 
disp{[’ as = ',num2str(aS)," So/Xt for root transfer matrix']) 


disp([’ Ki] “mun str (KI]),' 1b/in deflectiom at the Tone 41 


plotstat(2,rf,state,omega(1)) 

EXps[I1; 058050]; 
state=mykloop(tip,omega(1),rf,muf, EI£,Ohm) ; 
bTh=state(1,1); bX=state(2,1); bM=state(3,1); 
plotstat(3,rf,state,omega(1)) 


ESEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 
load h66cuff.mat -ascii 


nel=length(h66blade(:,1)); 

rcsfliplrthé6cuff(:,l)'*raduusie $ in 
wcsEPTDIPOR66cuff(:,2)9)5 z Thin? 
razrb(17:length(rb)); 

wazwb(17:length(wb)); 


rca=[rc,ra]; 

$ in 
wca=[wc, wa]; 

3 Than? 


S=0;I=0;M=0; 

for n=2:length(rca) 
dr=rca (n) -rca (n-1); 
wdr=mean (wca (n-1:n)); 
rcent=mean (rca (n-1:n)); 
S=S+rcent*dr*wdr; 
I=I+rcent"*2*dr*wdr; 

=M+dr *wdr ; 
end 


pgb-(e152*5M^2-2*e1*M*S*:8^2)7(5e1^2*M 59e 1* SED) > 
MS-IS^2ZSIT*M)/(-e1^2*M*2*e1*S-1); 
B= (el I1)/ (-el*’Mrs)e 


disp([’ Mb = ‘,num2str(Mb),’ slugs for blade mass at position R']) 


disp([' Ms = ',num2str (Ms),' slugs for mass at snubber position']) 
disp([' R s ',num2str(R),' inches radial location for mass Mb']) 
drsp(' ^) 


EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 


function state=mykloop(tip,omega,r,mu, EI, Ohm) 
% next step subfunction in Myklestad method 
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nbe-length (r); 


12[diff(F)]:120d 2 (Length) we $ delta-r in 
V=IW/ EI; 
% 1/tin*lb) 
Uz(1.^2) WN2*EI)S 
1/15 
Gz(l.^3 9 70S * ET) 
dnd 
Tsfiiplir(cumsum(flipir(mu.*r)))*Ohm^2; $ lb 


$ tip correction to shear: 
tip(4)=tip (4) +(mu (nbe-1)*1 (nbe-1) * (omega*2+Ohm*2)*tip(2))*.47; 


state=[zeros(4,nbe-1),tip]; % initialize state 


for n-(nbe-1):-1:1 


Th=state(1,n+1); $ theta 
X=state(2,n+1); % deflection 
M=state(3,n+1); $ moment 
S=state(4,n+1); % shear 


Thn = Th*(1+T(n)*U(n))-M*V(n)-S*U(n); 

Xn - X-Thn*l(n)-T(n)*Th*G(n)-M*U(n)-S*G(n); 
Mn = M+S*1 (n) -T (n) * (X-Xn) ; 
Sn = S+mu(n)*1(n)* (omega*2+Ohm*2) *Xn; 


$ [Thn; Xn;Mn; Sn] 
debug 


state(:,n)=[Thn;Xn;Mn;Sn]; 
end 
$state 


ESSELESELELELEEEEELEEEEEEEEEEEEEEESEESEEEEEEEEEEEEEEEEEE SES 
function [omegan, dynam, state] =natural(r,mu,EI,Ohm,articu,toler) 


9 


$ Finds natural frequency and dynamics matrix using 


Le 


% Myklestad method for beam bending 


% Robert L. King 
April 1998 


oe 


k=0; 
nbe=length(r); 


if nbe-=length (EI) 


error('The El vector is not the same size as the r vector.') 
elseif nbe--1length (mu) 
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error('The mass/length vector is not the same size as the r 
vector.) 
end 


if articu==1 t first guess at 1st bending freq - use eqn for 

uniform beam: 
omega_1=15.4*sqrt (mean (EI) / (mean (mu) * (r (l1length(r))-r(1))74)); 
% rad/sec 

else 
omega_1=3.52*sqrt (mean(El) / (mean (mu) * (r (length (r))-r(1))74)); 
% rad/sec 

end 

omega=[0.9 1.1]*omega_1; 


Qo ckckck ck ck kk k k kk k*k*k*WHTLE Role)’. 2.2.2.2. 2.2.2 2 2 2 2 2. 2 2 2 2 2 2 2.2. 
while (abs (omega (1) -omega (2) )>toler) 
k=k+1; 


$ Th tip - 1, all others zero 
tib-[(17 0; 070]; 
state-mykll(tip,omega(1),r,mu,EI,Ohm); 


% compute b components of dynamics matrix 
bTh=state(1,1); bX=state(2,1); bM=state(3,1); 
$b= (bTh, bX, bM] 
% debug 
%plotstate (2,r,state,omega(1)),pause 
% debug 


$ X tip - 1, all others zero 
tipzporei; 0; 015 
statezmykll(tip,omega(1),r,mu,EI,Ohm); 


% compute a components of dynamics matrix 
aTh=state(1,1); aX=state(2,1); aM=state(3,1); 
$a=[aTh, aX, aM] 
% debug 
$plotstate(3,r, state, omega (1) ) pause 
$ debug 


$ compute determinant 
Jdiguartlicus- 


dynam- [bM,aM;bX,aX]; $ articulated blade 
else 
dynam= [bTh,aTh;bX,aX]; % cantilevered blade 


end % if statement 


dte 
omega-fliplr (omega); 
err=det (dynam) ; 
else 
err= [det (dynam) ,err] ; 
% use linear interpolation for next guess at omega 
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omega- [max(0,omega(1)-err(1)*((omega(2)-omega(1))/(err(2)- 


err(1)))), omega]; 
%[spline(err,omega,0), omega]: % optional spline 
method 
end 


end % while loop 


XARXA KK A TO I, 5 LOOP***-**** xxx kk kk kk k kk 


disp([' ',;num2Str(Kk)," iterations in natural.m.']) 
$ debug 
omegan-omega(1); 


z4 debugging 
figure (1) plot(omega(2:k+1), err,omega(2:k+1),err,'o') 


xlabel('omega'),ylabel('determinant'),grid,title('Convergence Plot') 
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APPENDIX D - MAPLE WORKSHEET FOR COLEMAN MODEL 


EQUATIONS OF MOTION FOR GROUND RESONANCE CONSIDERING ONLY INPLANE 
DEGREES OF FREEDOM 


se 


DEFINE COORDINATE TRANSFORMATIONS 


so 


Blade undeformed axis -- Hub : 
> restart: 
> with(linalg): 
Warning, new definition for norm 
Warning, new definition for trace 
> psi:zOmega"ttPhi[k]; 
psi :- Omega t + Phi[k] 


> Tl:=alpha->matrix(3,3,[1,0,0,0,cos(alpha), sin(alpha),0,- 
sin(alpha),cos(alpha)]); 


T1 : alpha -> matrix(3, 3, 
[1, 0, 0, 0, cos(alpha), sin(alpha), 0, -sin(alpha), cos (alpha) ]) 


> T2:=alpha->matrix(3,3, [cos (alpha) ,0,- 
sin(alpha),0,1,0,sin(alpha),0,cos(alpha)]): 


12 := alpha => matrix(3, 3, 
[cos (alpha), 0, -sin(alpha), 0, 1, 0, sin(alpha), 0, cos(alpha)]) 


> T3:=alpha->matrix(3,3, [(cos(alpha),sin(alpha),0,- 
sin(alpha),cos(alpha),0,0,0,1]); 


T3 :- alpha -» matrix(3, 3, 
[cos (alpha), sin(alpha), 0, -sin(alpha), cos(alpha), 0, 0, 0, 1]) 
> diffl:=arg->map (diff,arg,t); 
diff1 := arg -> map(diff, arg, t) 


> Ml:<transpose(T3(psi)); 


[cos (%1) -sin(%1) 0] 
M := en cos ($1) ni 
i 0 0 i 

%1 := Omega t + Philk] 
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> M2:=transpose(T3(zeta[k] (t))); 


[cos (zeta [X] (t)) -sin(zeta[k](t)) 0] 
[ ] 
M2 := [sin(zeta[k](t)) cos (zeta[k] (t) ) 0] 
[ ] 
[ 0 0 1] 


Energy of rotor blades 


Kinetic energy of kth rotor blade (TBk) 
» rhoHI 1:=vector((u[1) (t), ul2) (6), 01): 
rhoHr-I = a eon 
> rhoBuH:zvector([e1,0,0]); 
rhoBuH := [el, 0, 0) 
> rhoBuH I:zmultiply(M1,rhoBuH) ; 
rhoBuH I :- [cos(Omega t 4 Phi[k]) el, sin(Omega t t Phi[k]) el, 0] 
» rhoPBd:-vector([R,0,0]); 
rhoPBd := [R, 0, 0] 
> rhoPBd_I:=multiply(M1,M2,rhoPBada); 
rhoPBd I :- [(cos($1) cos(zeta[k](t)) - sin($1) sin(zeta[k](t))) R, 
(sin(%1) cos(zetalk](t)) + cos($1) sin(zetalk] (t))) R, 0] 
$1 := Omega t + Philk] 
> rho:=map (simplify,matadd (rhoHI. I,matadd(rhoBuH I,rhoPBd I))): 
rho += [ull] (E) + cos($1) el + Reees(%1) cos(zetalk] (t)) 
- R sin($1) sin(zeta[k] (t)), u[2] (t) * sin($1) el 
+ R sin(%1) cos(zeta[k)(t)) # R a sin(zeta[k](t)), 0] 
$1 := Omega t + Phi[k] 


5 V:cdifflírho); 


[/d \ 
vV := [|-- u[1](t)| - sin(%2) Omega el 
[Ndt / 
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- R sin(%2) Omega cos(zeta[k](t)) - R cos(%2) 
- R cos(%2) Omega sin(zeta[k](t)) - R sin(%2) 
ja N 
, |-- u[2](t)| + cos(32) Omega el 
dt / 
+ R cos(%2) Omega cos(zeta[k] (t)) - R sin(%2) 
- R sin(%2) Omega sin(zeta[k](t)) 4 R cos(%2) 
] 
220] 
d 
$1 := -- zetalk] (t) 
dt 
+2 := Omega t + Phi[k] 
ENEI EK :—1/2*mb[Xk])*(V[1]^2:*V[2]^2); 
///da N 
TBk :- 1/2 mb[k] |||-- ut1](t)| - sin($2) Omega ei 
NN AGE / 
- R sin(%2) Omega cos(zeta[k](t)) - R cos(%2) 
- R cos(%2) Omega sin(zeta[k](t)) - R sin(%2) 
\2 Fa N 
| + ||-- uf2](t)] + cos(%2) Omega el 
/ \\dt / 
+ R cos($2) Omega cos(zeta[k](t)) - R sin($2) 
~ R sin(%2) Omega sin(zeta[k] (t)) + R cos(%2) 
NAN 
| 
/ 
d 
41 := -- zetalk) (t) 
dt 
$2 := Omega t + Philk] 
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sin(zeta[k] (t)) 


cos (zeta[k] (t)) 


sin(zeta[k] (t)) 


cos (zeta [k] (t)) 


sin(zeta[k] (t)) 


cos (zeta[k] (t) ) 


Sin(zeta[k] (t) ) 


cos (zeta[k] (t)) 


$1 


$1 


$1 


$1 


$1 


$1 


%1 


$1 


S UBk1:z1/72*Ke[k]*zeta[k](t)^2; #Linear Elastic Forces 


2 
UBk1 :- 1/2 Ke[k] zeta[k](t) 


» UBk2:-1/4*Kd[k]*zeta[k] (t)^4; *Duffing Elastic Forces 


4 
UBk2 := 1/4 Kd[k] zeta[k] (t) 


> UBk3:=1/4*Ks[k] *signum(zetalk] (t)-z)*(zetalk] (t) *2+z”2- 
2*zeta[k] (t) *z) *1/4*Ks[k] *signum(zeta[k](t)-«z)*(-zeta[Xk] (t)^2-z^2- 
2*zeta[k] (t)*z)+1/2*Ks[k]*zeta[k] (t)^2+1/2*Ks[k]*z^2; 


UBk3 :z 
1/4 
2 2 
Ks[k] signum(zeta[k](t) - z) (zeta[k](t) +z - 2 zetal[k](t) z) 
* 1/4 
2 2 
Ks[k] signum(zeta[k] (t) + Zz) (-zeta[k] (t) -Z - 2 zeta[k](t) z) 
2 2 
439172 Ks[k] zetalk](t) + 1/2 Ks[k] z 
> UBk:=UBk1+UBk2+UBk3; 
> UBk:=UBk1+UBk2; 
2 4 


UBk := 1/2 Kelk) zetalk](t) + 1/4 Kd[k] zeta[k] (t) 


em wm wm wm wm em em wm wm vm wm vm vm wm ge vm ge ge mm vm ge wm ge emm gem mm vm e em wm gn wm En ke m En En En 


DBk:=1/2*Czeta[k]*(dđiff(zeta[k](t),t))^2+Vzeta[k]*(điff(zeta[k](t),t))^ 
2>abs(ditf(zetalk])(t),E)): 


/a NO 
DBk := 1/2 Czeta[k] |-- zeta[k] (t) | 
\dt / 
/d A: 
* Vzeta[k] |-- zeta[k](t)| | -- zeta[k] (€) | 
Nat Mae | 


Energy of hub 
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Kinetic energy of hub (TH) / Potential energy of hub (UH) / Dissapative 
function of hub (DH) 


= TE:=1/2*M11]*(4xftf£(u114 (6) t))o2e1/2*MI2] EE 


/a \2 /d N2 
TF :- 1/2 M({1) |-- ull] (t)| + 1/2 M[2] |-- ul2] (t) | 
Nat / Nat 7 


AUNE =17/ 25 LI] *u UI CE] 2412K [2G UZ] (CT) 22. 


2 2 
UP := 1/2 Khi] u[il(t) + 1722502 NN t) 


> 
DE:-1/2*e6[1]*(diff(u[i](t),t))^241/2*c[2] (41 E) (6) E 0002412 v [1] 
(a(t), 6))°2*ahs(d1ff (ul1] (t). t)) 172511 (Gt (loo ^2*3 
EE EG E ENK 


/d NG /d NG 
DF :< 1/2 c[1] |-- u[1](t)| + 1/2 ¢[2] |-- ul2] (t) | 
\dt Z \dt / 
/d \2 | a 
4 1/2 vI1] |== ult) (e)| ) | == ey 
AdE AdE 
/d 120 Hd 
+ 1/2 v[2] |-- u[2](t)| I] -- u[2]í(t) | 
\dt / | dt 


Generalized forces on generalized displacements 


> F[1]:=0; 
e[1] 9 
> F[2]:=0; 
ERI. 
ER UI Ez 
ED: - ur] 
> F[4]:=u[2]; 
F[4] := u[2] 
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> KSI- =u; 


RIS wp) 


Derivation of equations of motion using Lagrange's equation 


pORR:=[u11]ME) 121 (ES 


> 
> DOFB:-2[zeta[1] (t),zeta[2] (t),zeta[3] (t) ]: 
» DOF:-2[op(DOFF),op(DOFB)]: 
> Q@DOF:=diff1 (DOF): 
> dADOF:=diff1 (ADOF) : 
> setA:={}:setB:={}:setC:={}: 
> setD:={}:setE:={}:setF:={}: 
> DOFq:=[] :dDOFq:=[] :ddDOFq:=[]: 
> for i from 1 to vectdim(DOF) do 
> BoR@a:=-[op (DORF), a[i]]: 
> dDOFa:- [op(dpDOFa) , da[i]] : 
> ddDOFq:-2[op(ddDOFq),ddq[1ill: 
> setA:=setA union (ddDOF[i]-2ddDOFq[i]): 
> setB:=setB union {dDOF[i]=dDOFq[i]}: 
> setC:=setC union {DOF[i]=DOFq[i]}: 
> setD:=setD union {ddDOFq[i]=ddDOF[i]}: 
> setE:=setE union {dDOFq[i]=dDOF[i]J}: 
> setF:=setF union {DOFq[i]=DOF[i]}: 
> od: 
> setl:=setA union setB union setC; 
2 
d d 
secl ---— ult) ddatil, -— ulj (3) = da] rule (E) apri 
2 dt 
de 
2 
d d 
EE, ddgi2], -- u[21(t) - dq[2], u[2] (5) - a[215 
2 dt 
dt 
2 
d d 
--- zeta[1](t) » ddq[3], -- zeta[1](t) < dg[3], zetal1](t) < gl3], 
2 dt 
dt 
2 
d d 
--- zeta[2](t) = ddq[4], -- zeta[2](t) = dq[4], zeta[2](t) = q[4], 
2 dt 
dt 
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2 

d d 

--- zeta[3](t) s ddq[5], -- zetal3] (t) - da[(5], zeta[(3](t) = ql5] 
2 dt 

dt 


) 


> set2:=setD union setE union setF; 


2 
d d 
Seca se tadg[1] =, -= 14106) dqd[i] = YAN 242 Je= uie) 
2 dt 
dt 
2 
d d 
dae = s2- w2 el; CEl2] ee 
2 dt 
dt 
2 
d d 
ddq[3] zs --- zeta[1](t), dq[3] » -- zeta[1](t), a(3] = zeta[1] (t), 
2 de 
GE 
2 
d d 
ddg[4] = --- zeta[2] (t), da[4] » -- zeta[2] (t), ql4] = zetal2] (t), 
2 dt 
dt 
D 
d d 
dda[5] = === zetal3](t), da[(5] - -- zeta[3](t), a(5] - zeta[3](t) 
2 dt 
AE 
) 
E: 
U:=UF: 
D1:=DF: 


VV VVV VV VV VVV 


for 1 £rom l to vectdim(DOFB) do 
T:=T+subs(k=i,TBk): 
U:=U+subs (k=i, UBk) : 
D1:=D1+subs (k=i,DBk): 
od: 
Temp :=subs(set1,T): 
for i from 1 to vectdim(DOF) do 
temp1:2diff(Temp,dDOFq[i]): 
temp2:=subs (set2,templ): 
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> temp3 :=diff (temp2,t): 

> L1:2subs(setl,temp3): 

> L2:=diff(Temp,DOFq[i]): 

> L3:<diff(subs(set1,U),DOFgli]): 

> L4:=diff (subs (set1,D1), dDOFq[i]): 

> EOM[i] :=simplify(L1-L2+L3+L4-F[i]): 
> od: 


> setS:={signum(1,q[3]-z)=0,signum(1,q[4]-z)=0,signum(1,q[5]- 
z)=0,signum(1,q[3]+z)=0,signum(1,q[4) +z) =0,signum(1,q[5]+z)=0,abs(1,dq[ 
3] }=0,abs(1,dq[4])=0,abs(1,dq[5])=0,abs(1,dq[1])=0,abs(1,dq[2] )=0}: 

> A:zmatrix(vectdim(DOF) ,vectdim(DOF)); 


Afi= array(l NE MA). SPE] 


Bor i from 1 to vectdim(DOF) do 
for j from 1 to vectdim(DOF) do 
A[i,j]:=coeff(EOM[i],ddDOFq[j]): 
A[1,3]:=subs(setS,A[1,3]): 
od: 
od: 
Ax2dot:=multiply (A, ddDOFq) : 
f:zarray(l..vectdim(DOF)); 


V V V V OVO VOV V 


Emo sarravii".. 5, .1]) 


for i from 1 to vectdim(DOF) do 
f [i] :=-simplify(EOM[i]-Ax2dot[i]): 
£{i]:=subs(setS,f[i]): 
od: 
zidot:z(|rexi:: | |: 
for i from 1 to vectdimi(DOE) do xlIdot:s[op(xldot)?x[i]l] od: 
for i from vectdim(DOF)+1 to 2*vectdim(DOF) do x1:=[op(x1),x[i]] od: 
setX:=(): 
for i from 1 to vectdim(DOF) do 
setX:=setX union (dDOFq[i]=x1dot[1]): 
setX:=setX union (DOFq[1]=x1[1]): 
od: 
interface(labelling-false); 
Al:-subs(setX ,op(A)); 


V VV V V V V V NV NV V N N N 


> 
ren 


MAA O 95 
-mb[1] R cos($1) sin(x[8]) - mb[1] R sin($1) cos(x[8]) , 
-mb[2] R cos($2) sin(x[9]) - mb[2] R sin(32) cos(x[9]) | 
-mb[3] R cos($3) sin(x[10]) - mb[3] R sin($3) cos(x[10])] 
[0 , MI2] + mb[1] + mb[2] + mb[3] , 


-mb[1] R sin(%1) sin(x[8]) + mb[1] R cos(%1) cos(x[8]) , 
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-mb[2] R sin(%2) sin(x[9]) + mb[2] R cos(%2) cos(x[9]) , 
-mb[3] R sin($3) sin(x[10]) + mb[3] R cos(%3) cos(x[10])] 


[ 
[-mb[1] R cos($1) sin(x[8]) - mb[1] R sin(%1) cos(x[8]) , 


2 
-mb[1] R sin($1) sin(x[8]) + mb[1] R cos($1) cos(x[8]) , mb[1] R 
AD O) 


[ 
[-mb[2] R cos(%2) sin(x[9]) - mb[2] R sin($2) cos(x[9]) , 


-mb[2] R sin($2) sin(x[9]) + mb[2] R cos(%2) cos(x[9]) , 0, 


2 ] 
mb[2] R , 0] 


[ 
[-mb[3] R cos($3) sin(x[10]) - mb[3] R sin($3) cos(x[10]) , 


-mb[3] R sin(%3) sin(x[10]) + mb[3] R cos(%3) cos(x[10]) , 0, 0 


2] 
pump a 


$1 := Omega t + Phi[1] 
$2 := Omega t + Phi[2] 
$3 := Omega t + Phi[3] 
> £1:=subs(setX ,op(£)); 


[ 
aT: [(SK[Z] x[6] = €[1] x[1] - vlil x[1] | x[1] | 


2 2 
+ mb[1] cos(%2) Omega el + mb[1] R cos(%2) Omega cos(x[8]) 


~ 2 mb[1] R sin(%2) Omega sin(x[8]) x[3] 


2 
+ mb[1] R cos(%2) cos(x[8]) x[3] 


2 
- mb[1] R sin(%2) Omega sin(x[8]) 


+ 2 mb[1] R cos(%2) Omega cos(x[8]) x[3] 


2 2 
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mb[1] R sin($2) sin(x[8]) x[3] + mb[2] cos(%1) Omega el 


2 
mb[2] R cos($1) Omega cos(x[9]) 


2 mb[2] R sin($1) Omega sin(x[9]) x[4] 


2 
mbiž] R cos(%1) cos(xi9]) x[4] 


2 
mb[2] R sin($1) Omega sin(x[9]) 


2 mb[2] R cos($1) Omega cos(x[9]) x[4] 


2 2 
mb[2] R sin($1) sin(x[9]) x[4] + mb[3] cos($3) Omega el 


2 
mb[3] R cos($3) Omega cos(x[10]) 


2 mb[3] R sin(%3) Omega sin(x[10]) x[5] 


2 
mb[3] R cos(%3) cos(x[10]) x[5] 


2 
mb[3] R sin($3) Omega sin(x[10]) 


2 mb[3] R cos($3) Omega cos(x[10]) x[5] 


2 
MS A szdatóz) Sin (2:1 191) x[(5] . -e[(2] x[2] 2 KIZ] #[7] 


2 
v[2] x[2] | x[2] | + mbi2] R cos(%1) Omega sin(x[9]) 


2 mb[2] R sin(%1) Omega cos(x[9]) x[4] 


2 2 
mb[2] R cos(%1) sin(x[9]) x[4] + mb[3] sin(%3) Omega el 


2 
mb[3] R sin(%3) Omega  cos(x[10]) 


2 mb[3] R cos($3) Omega sin(x[10]) x[5] 


2 
mols] R sin($3) cos(x[10]) x[5] 


2 
mb[3] R cos($3) Omega sin(x[10]) 


2 mb[3] R sin($3) Omega cos(x[10]) x[5] 
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2 


+ mb[3] R cos/(83) sIn(#[10]) x[5] 
2 
+ mb[1] R sin($2) Omega cos(x[8]) 


+ 2 mb[1] R cos(%2) Omega sin(x[8]) 
2 
+ mb[1] R sin($2) cos(x[8]) x[3] 
2 
+ mb[1] R cos(%2) Omega sin(x[8]) 


+ 2 mb[1] R sin($2) Omega cos(x[8]) 
2 
+ mb[1] R cos(%2) sin(x[8]) x[3] + 
2 
+ mb[2] R sin(%1) Omega cos(x[9]) 
2 
+ mb[2] R sin(%1) cos(x[9]) x[4] 


+ 2 mb[2] R cos(%1) Omega sin(x[9]) 


vzeta[1] x[3] | x[3] 
3 
- Ke[1] x[8] - Kd[1] x[8] 


- Czeta[2] x[4] - Ke[2] x[9] 
2 
- mb[2] Omega el R sin(x[9]), 


3 
=shels) x[10)"— Kalo] ex (10) 
2 ] 


- mb[3] Omega el R sin(x[10])] 


$1 :- Omega t 4 Phi[2] 
$2 := Omega t + Phi[1] 
$3 := Omega t + Phil3] 


> realib(fortran): 
> B:=augment(A1,£1): 
> fortran(B, optimized) ; 
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+ mb[1] 


| - mb[1] Omega 


- Czetal1] x[3], 


- 2 Vzeta[2] x[4] 


-2 Vzeta[3] x[5] 


2 
sin($2) Omega el 


[343 


XIS] 


2 


mb[2] sin($1) Omega el 


ac [415 


2 

el R sin(x[8]) + u[1] 
8 
-Kd[2] x[9] 


| x[4] | t u[2] 


| x[5] | + u[3] 


- Czeta[3] x[5] 








APPENDIX E - EQUATIONS OF MOTION FOR COLEMAN MODEL 


mb[3]*R*sin(Omega*t4Phi[3])*sin(q[5]) *dq[5]^2- 
mb[3]*R*sin(Omega*t+Phi[3])*cos(q[5])*ddq[5]+mb[2]*R*sin(Omega*t+Phi[2])*Ome 
ga^2*sin(q[4])- 
2*mb[2]*R*cos(Omega*t+Phi[2])*Omega*cos(g[4])*dg[4]+mb[2]*R*sin(Omega*t+Phi 
[2])*sin(q[4]) *dq[4]^2-mb[2] *R*sin(Omega*t--Phi[2])*cos(q[4  *ddq[4]- 
mb[3]*cos(Omega*t4Phi[3])*Omega^2*e1- 
mb[3]*R*cos(Omega*t+Phi[3])*Omega*2*cos(q[5])+2*mb[3]*R*sin(Omega*t+Phi[3]) 
*Omega*sin(q[5])*dq[5]-mb[3]*R*cos(Omega*t+Phi[3])*cos(q[5])*dq[5]*2- 
mb[3]*R*cos(Omega*t+Phi[3])*sin(q[5])*ddq[5]+mb[3]*R*sin(Omega*t+Phi[3])*Ome 
ga^2*sin(q[5])- 
mb[1]*R*cos(Omega*t+Phi[1])*Omega’2*cos(q[3])+2*mb[1]*R*sin(Omega*t+Phi[1]) 
*Omega*sin(q[3]) *dq[3]-mb[1]*R*cos(Omega*tPhi[1])*cos(q[3])*dq[3]^2- 
mb[1]*R*cos(Omega*t+Phi[1])*sin(q[3])*ddq[3]+mb[1]*R*sin(Omega*t+Phi[1])*Ome 
ga^2*sin(q[3]- 
2*mb[1]*R*cos(Omega*t+Phi[1])*Omega*cos(q[3])*dq[3]+mb[1]*R*sin(Omega*t+Phi 
[1])*sin(q[3])*dq[3]42-mb[1]*R*sin(Omega*t+Phi[1])*cos(q[3])*ddq[3]- 
mb[2]*cos(Omega*t+Phi[2])*Omega’2*el - 

mb[2] *R *cos(Omega*t+Phi[2]) *Omega'2*cos(g[4])+2*mb[2]*R*sin(Omega*t+Phi[2]) 
*Omega*sin(q[4])*dq[4]-mb[2]*R*cos(Omega*t-Phi[2])*cos(q[4])* dq[4]^2- 
mb[2]*R*cos(Omega*t+Phi[2])*sin(q[4])*ddq[4]- 
mb[1]*cos(Omega*t+Phi[1])*Omega‘’2*e1- 
2*mb[3]*R*cos(Omega*t+Phi[3])*Omega*cos(q[5])*dq[5]+K[1]*q[1]+v[1]*dq[1]*abs( 
dg[1])+M[1]*ddgq[1]+1/2*v[1]*dq[1]*2*abs(1,dq[1])+c[1]*dq[1]+mb[1]*ddq[1]+mb[2]* 
ddq[1]+mb[3]*ddq[1]=0 


mb[1]*ddq[2]+mb[3]*ddq[2]-mb[3]*R*cos(Omega*t+Phi[3])*FOmega”2*sin(q[5])- 
2*mb[3]*R*sin(Omega*t+Phi[3])*Omega*cos(q[5])*dq[5]- 
mb[3]*R*cos(Omega*t+Phi[3])*sin(q[5])*dq[5]42+mb[3]*R*cos(Omega*t+Phi[3])*cos 
(q([5])*ddq[5]-mb[3 ]*R*sin(Omega*t+Phi[3])*Omega*2*cos(q[5])- 
2*mb[3]*R*cos(Omega*t+Phi[3])*Omega*sin(q[5])*dq[5]- 
mb[3]*R*sin(Omega*t+Phi[3])*sin(q[5])*ddq[5]- 
2*mb[2]*R*sin(Omega*t+Phi[2])* Omega*cos(q[4])*dq[4]- 
mb[3]*R*sin(Omega*t+Phi[3])*cos(q[5])*dq[5]42- 
mb[1]*sin(Omega*t«Phi[1])*Omega^2*e1- 
mb[1]*R*sin(Omega*t*Phi[1])*Omega^2*cos(q[3])- 
2*mb[1]*R*cos(Omega*t4-Phi[1])*Omega*sin(q[3])*dq[3]- 
mb[1]*R*sin(Omega*tPhi[1])*cos(q[3])*dq[3]^2- 
mb[1]*R*sin(Omega*t«Phi[1])*sin(q[3])*ddq[3]- 
mb[1]*R*cos(Omega*t-4-Phi[1])*Omega^2*sin(q[3])- 
2*mb[1]*R*sin(Omega*tPhi[1])*Omega*cos(q[3])*dq[3]- 
2*mb[2]*R*cos(Omega*t+Phi[2])*Omega*sin(q[4])*dq[4]- 
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mb[2]*R*sin(Omega*t+Phi[2])*cos(q[4])*dq[4]*2- 
mb[2]*R*sin(Omega*t+Phi[2])*sin(q[4])*ddq[4]- 
mb[2]*R*cos(Omega*t+Phi[2])*Omega'2*sin(g[4])- 
mb[1]*R*cos(Omega*t+Phi[1])*sin(q[3])*dq[3]42+mb[1]*R*cos(Omega*t+Phi[1])*cos 
(q[3)*ddq[3]-mb[2]*sin(Omega*t4-Phi[2]) *Omega^2*e1- z 
mb[2]*R*sin(Omega*t4-Phi[2])*Omega^2*cos(q[4])- 

mb[3]*sin(Omega*t4-Phi[3]) *Omega^2*e1- 
mb[2]*R*cos(Omega*t+Phi[2])*sin(q[4])*dq[4]*2+mb[2]*R *cos(Omega*t+Phi[2])*cos 
(q[4])*ddq[4]+M[2]*ddq[2]+mb[2] *ddq[2]+1/2*v[2]*dq[2]42*abs(1,dq[2])+K[2]*q[2]+ 
c[2]*dg[2]+v[2]*dg[2]*abs(dg[2])=0 


2*Vzeta[1]*dq[3]*abs(dq[3])- 
u[{1]+mb[1]*ddg[2]*R*cos(Omega*t+Phi[1])*cos(q[3])+mb[1]*R%2*ddq[3]- 
mb[1]*ddg[1]*R*cos(Omega*t+Phi[1])*sin(g[3])- 
mb[1]*ddq[1]*R*sin(Omega*t+Phi[1])*cos(q[3])- 
mb[1]*ddq[2]*R*sin(Omega*t+Phi[1])*sin(q[3])+mb[1]*Omega’2*e1*R*sin(q[3])+Vze 
ta[1]*dq[3]42*abs(1,dq[3])+Ke[1]*q[3]+Kd[1]*q[3]*3+Czeta[1]*dq[3]=0 


mb[2]*Omega^2*e1*R*sin(q[4])- 
mb[2]*ddg[1]*R*sin(Omega*t+Phi[2])*cos(q[4])+mb[2]*ddq[2]*R*cos(Omega*t+Phi[2 
)*cos(q[4])2-mb[2] *R^2*ddq[4]- 

mb[2]*ddq[2]*R *sin(Omega*t4-Phi[2 ])*sin(q[4])- Vzeta[2] *dq[4]^2*abs(1,dq[4])- 
u[2]+Ke[2]*q[4]+Kd[2] *q[4]*3+Czeta[2]*dq[4]+2* Vzeta[2]*dq[4]*abs(dq[4])- 
mb[2]*ddq[1]*R*cos(Omega*t+Phi[2])*sin(q[4])=0 


-u[3]-mb[3]*ddq[1]*R*cos(Omega*t--Phi[3])*sin(q[5])- 
mb[3]*ddg[1]*R*sin(Omega*t+Phi[3])*cos(g[5])+mb[3]*Omega'2*e1*R*sin(g[5])+mb 
[3]*RA2*ddg[5]- 
mb[3]*ddq[2]*R*sin(Omega*t+Phi[3])*sin(q[5])+mb[3]*ddq[2]*R*cos(Omega*t+Phi[3 
])*cos(q[5])+2* Vzeta[3]*dq[5]*abs(dq[5])+Vzeta[3]*dq[5]42*abs(1,dq[5])+Ke[3]*q[5]+ 
Kd[3]*q[5]^34-Czeta[3] *dq[5]20 
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APPENDIX F - SIMPLIFIED ROTOR/FUSELAGE EQUATIONS OF MOTION 


Equations with viscous and hydraulic damping in the hub 
degrees of freedom, x and y. [Ref. 19] 


EE 

— mb, * R* (C, * sin, - £,) - (Q-- £,)? *cos(V +2) 
t mb, * R* (C, *sin(V, +,)+ (A +C,)? *cos(YP, +£,) 
mb, * R* (C, *sin(E, + $3) + (Q+ £,)? *cos(P, 9 £,) 


M,*j4C,* y*V, * y*ly| - K,*y 

— mb, * R* (C, "sin(E, - £,) - (Q-- £,)? *sin(E, ő) 
t mb, * R* (C, *sin(, - £,) -(Q- £,)? *sin(VP,  £,) 
mb, RAC Fsin(W, E22) (0405) #sin(F, EE 


mb, *R? s +C *C, +Ke, st, * mb, O" *e*R*sin(6,) 
= mb, * i * R*sin(Y, FO) mb *y*R*cos(P, +61) 


mb, RO sle TC kd, t Ke, #C,t mb, O" "e"R"sin(4,) 
= mb, * X* R*sin( WP, + (,) = mb, = y R®cos(h, +¢,) 


mb, *R? M no SE +Ke, * €, + mb, * 2° *e*R *sin(/,) 
= mb, 7X Rost, +63) — mb, t Ý * R *cos(F, GO 
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APPENDIX G - 5 BLADED EQUATIONS OF MOTION WITH POLYNOMIAL 
SNUBBER FUNCTIONS 


2 
mb[4] R sin(Omega t + Phi[4]) Omega sin(q[6]) 


- 2 mb[4] R cos(Omega t + Phi[4]) Omega cos(q[6]) dq[6] 


2 
4 mb[4] R sin(Omega t t Phi[4]) sin(gl6]) dal[6] 


- mb[4] R sin(Omega t t Phi[4]) cos(g[l6]) ddal[6] 


2 
- mb[5] cos(Omega t + Phi[5]) Omega el 


2 
- mb[5] R cos(Omega t + Phi[5]) Omega cos(q[7]) 


+ 2 mb[5] R sin(Omega t + Phi[5]) Omega sin(g[7]) dq(7] 
- 2 mb[5] R cos(Omega t + Phi[5]) Omega cos(q[7]) dq[7] 


2 
- mb[2] R cos(Omega t + Phi[2]) cos(g[4]) daąal4] 


- mb[2] R cos(Omega t + Phi[2]) sin(q[4]) ddq[4] 


2 
+ mb[2] R sin(Omega t + Phi[2]) Omega sin(q[4]) 


- 2 mb[2] R cos(Omega t + Phi[2]) Omega cos(a[4]) da[4] 


2 
- mb[2] R sin(Omega t + Phi[2]) sin(q[4]) dq(4)] 


- mb[2] R sin(Omega t + Phi[2]) cos(q[4]) ddq[4] 


E 
- mb[3] R cos(Omega t + Phi[3]) Omega cos(q[5]) 


+ 2 mb[3] R sin(Omega t + Phi[3]) Omega sin(q[5]) dq[5] 


2 
- mb[3] R cos(Omega t 4 Phi[3]) cos(g[5]) dal[5] 


- mb[3] R cos(Omega t + Phi[3]) sin(g[5]) ddq[5] 


2 
+ mb[3] R sin(Omega t + Phi[3]) Omega sin(q[5]) 


2 mb[3] R cos(Omega t + Phi[3]) Omega cos(q[5]) dq[5] 
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mb[1] cos(Omega t + 


TO (OR 


2 mb[1] 


mb[1] R 


mb[1] R 


mb[1] R 


2 mo[1] 


mb[1] R 


mb[1] R 


cos (Omega t 


R sin(Omega 


cos (Omega t 


cos (Omega t 


sin(Omega t 


R cos (Omega 


sin(Omega t 


sin(Omega t 


mb[2] cos(Omega t + 


mb[2] R 


cos (Omega t 


mb[3] cos(Omega t + 


2 mb[2] 


mb[3] R 


mb[5] R 


mb[5] R 


mb[5] R 
mb[5] R 


mb[3] R 


R sin(Omega 


sin(Omega t 


cos (Omega t 


sin(Omega t 


sin(Omega t 
sin(Omega t 


sin(Omega t 


mb[4] cos(Omega t + 


mb[4] R 


cos (Omega t 


2 
Phi[1]) Omega el 


D 
* Phi[1]) Omega cos(q[3]) 


t * Phi[1]) Omega sin(q[3]) dq[3] 


2 
+ Phi[1]) cos(al31) da[3] 


4 Phi[1]) sin(a[3]) dda[3] 


2 
+ Phi[1]) Omega sin(a[3]) 


t + Phi[1]) Omega cos(q[3]) dq[3] 


2 
Eller lee Ed 


+ 


+ 


Phi[1]) cos(ql3]) ddq[3] 


2 
Phi[2]) Omega el 


2 
+ Phi[2]) Omega cos(q[4]) 


2 
Phi[3]) Omega el 


t + Phi[2]) Omega sin(q[4]) dq[4] 


2 
+ Phi[3)) sin(q[5]) dq[5] 


S PhI[5]) /sinigi7]) dgadi7) 


2 
+ Phi[5]) Omega sin(q[7]) 


2 
-4Ph1(51) sin(a171) 4317) 


4: Phi[5]) cos(a[7]) dda[7] 
3: Phi[3]) cos(q[5]) dda[5] 


2 
Phi[4]) Omega el 


2 
+ Phi[4]) Omega  cos(q[6]) 
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x [2] 


A 


2 mb[4] R sin(Omega t * Phi[4]) Omega sin(qd[6]) dq[6] 


2 
mb[4] R cos(Omega t + Phi[4]) cos(q[6]) dq[{6] - u[1] 


mb[5] ddq(1] * mb[4] ddq[1] 

mb[2] ddg[(1] + M[1] ddq[1] + mb[1] ddq[1] 

nid dali ara ECC CPI tt 
mb[4] R cos(Omega t + Phi[4]) sin(q[6]) ddq[6] 


2 
mb[5] R cos(Omega t + Phi[5]) cos(a[7]) dg[7] = 0 


datz al da 210] 


2 mb[3] R cos(Omega t + Phi[3]) Omega sin(q[5]) dq[5] 


2 
mb[3] R sin(Omega t + Phi[3]) cos(q[5]) dq[5] 


mb[3] R sin(Omega t + Phi[3]) sin(q[5]) ddq[5] 


2 
mb[3] R cos(Omega t + Phi[3]) Omega sin(q[5]) 


2 mb[5] R cos(Omega t + Phi[5]) Omega sin(q[7]) dq(7] 


2 
mb[5] R sin(Omega t + Phi[5]) cos(q[7]) da[7] 


mb[5] R sin(Omega t + Phi[5]) sin(a[7]) ddq[7] 


2 
mb[5] R cos(Omega t + Phi[5]) Omega sin(g[7]) 


2 mb[5] R sin(Omega t + Phi[5]) Omega cos(g[7]) dg[7] 


2 
mb[5] R cos(Omega t + Phi[5]) sin(q[7]) dq[7] 


mb[5] R cos(Omega t + Phi[5]) cos(q[7]) ddq[7] 
2 mb[3] R sin(Omega t * Phi[3]) Omega cos(q[5]) dq[(5] 


2 
mb[3] R cos(Omega t 4 Phi[3]) sin(g[5]) da[5] 


mb[3] R cos(Omega t + Phi[3]) cos(q[5]) ddq[5] 
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mb[4] sin(Omega t + 


mb[4] R 


2 mb[4] 


mb[4] R 


mb[4] R 


mb[4] R 


2 mb[4] 


mb[4] R 


mb[4] R 


mb[5] 


mb[1] 


mb[ TIME 


2 mb] 


mort] R 


mb[1] R 


mb[1] R 


2 mb[1] 


mb[1] R 


mb[1] R 


mb[2] 


mb[2] R sin(Omega t + Phi[2]) Omega 


sin(Omega t 


R cos (Omega 


sin(Omega t 


sin(Omega t 


cos (Omega t 


R sin(Omega 


cos (Omega t 


cos (Omega t 


sin(Omega t + 


sin(Omega t + 


sin(Omega t 


R cos (Omega 


sin(Omega t 


sin(Omega t 


cos (Omega t 


R sin(Omega 


cos (Omega t 


cos (Omega t 


sin(Omega t + 


2 
Phi[4]) Omega el 


2 
* Phi[4]) Omega cos(q[6]) 
t t Phi[4]) Omega sin(q[6]) dq[6] 


2 
+ Phi[4]) cos(q[6]) dq[6] 


+ Phi[4]) sin(q[6]) ddgq[6] 


2 
+ Phi[4]) Omega sin(q[6]) 
t + Phi[4]) Omega cos(q[6]) dq[6] 


2 
Phi[4]) sin(q[61) dq[6) 


+ 


+ 


Phi [4]) cos(q[6]) ddg[6] 


2 
Phi[5]) Omega el 


2 
Phi[1]) Omega el 


2 
+ Phi[1]) Omega cos(q[3]) 
t + Phi[1]) Omega sin(q[3]) dq[3] 


2 
4 Phi[1]) cos(a[(3]) dg[3) 


*wBhirl]) sin(q[3]) deég[3] 


2 
+ Phi[1]) Omega sin(q[3]) 
t + Phi(1]) Omega cos(q[3]) da[3] 


- 2 
it Phi[1]) sin(g[3]) agqs] 


+ Phi[1]) cos(q[3]) ddal3] 


2 
Phi[2]) Omega el 


2 
cos (q[4]) 
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2 mb[2] R cos(Omega t 4 Phi[2]) Omega sin(q[4]) da[4] 


2 
mb[2] R sin(Omega t + Phi[2]) cos(q[4]) da[4] 


mb[2] R sin(Omega t + Phi[2]) sin(q[4]) ddq[4] 


2 
mb[2] R cos(Omega t + Phi[2]) Omega sin(q[4]) 


2 mb[2] R sin(Omega t + Phi[2]) Omega cos(q[4]) dq[4] 


2 
mb[2] R cos(Omega t + Phi[2]) sin(q[4]) dq[4] 


mb[2] R cos(Omega t + Phi[2]) cos(q[4]) ddq[4] 


2 
mb[3] sin(Omega t + Phi[3]) Omega el 


2 
mb[5] R sin(Omega t + Phi[5]) Omega cos(q[7]) + mb[2] ddq[2] 


M[2] ddq[2] - u[2] + K[2] ql2] 


2 
c[2] dg[2] - mb[3] R sin(Omega t + Phi[3]) Omega cos(q[5]) 


mb[1] ddq[2] + mb[5] ddq[2] + mb[3] ddq[2] + mb[4] ddq[2] = 0 


- mb[1] ddg[2] R sin(Omega t + Phi[1]) sin(q[3]) 


+ 


mb[1] ddq[2] R cos(Omega t + Phi[1]) cos(q[3]) 
mb[1] ddq[1] R cos(Omega t - Phi[1]) sin(q[3]) 
mb[1] ddq[1] R sin(Omega t + Phi[1]) cos(q[3]) 


2 2 
mb[1] Omega el R sin(q[3]) + mb[1] R ddq[3] + q[3] Kpoly[5] 


de [31-^6polyt5] 


3 
a[3] Epoly[2] | abi] TEIS kost (AM al 
3 
dal3] epeolv[2] anl l m aa ere Eé E | ahil | 
2 4 5 


Eeer EE, ee i qs] S23 qls Ne ov m 
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+ 


zc 


2 


mb[2] Omega 


mb[2] 
mb[2] 
mb[2] 


mo 12] 


2 q[4] 


dq[4] Cpoly[5] + mb[2] R ddg[4] - u[4] + dq[4] Cpoly[1] q[4] 


a 
Kpoiky[3] - up3)---0 


el R sin(q[4]) 
ddq[1] R cos(Omega t + Phi[2]) sin(q[4]) 


ddq[1] R sin(Omega t 


4 


Phi[2]) cos(a[4]) 


ddq[2] R sin(Omega t 


4 


Phi[2]) sin(q[4]) 


ddq[2] R cos(Omega t 


+ 


Phi[2]) cos(a[4]) 


3 5 
Kpoly[3] + 3 g[4] Kpoly[1] + g[4] Kpol 


2 


3 
q[4] Kpoly[2] | al4] | + g[4] Kpoly[4] | gq[4] 
da[4] Cpoly[4] | a[4] | 
3 2 
dq[4] Cpoly[2] | a[4] | + dql4] Cpoly[3] al4] 
2 


mb[3] R ddq[5] 


+ 


mb[3] 
mb[3] 
mb[3] 


mb[3] 


mb[3] 


dq[5] 


aq] 


ddq[2] R sin(Omega t * Phi[3]) sin(q[5]) 
ddq[2] R cos(Omega t + Phi[3]) cos(g[5]) 
ddq[1] R cos(Omega t + Phi[3]) sin(q[5]) 
ddq[1] R sin(Omega t + Phi[3]) cos(q[5]) 


2 
Omega el R sin(g(5]) 


4 E 
Cpoly[i] al5] 9 = wisi ds (Siero ly 210] 


Cpelyi4) | qi5] (+ a(S] Keelyi4) (eis) 


3 2 


q[5] Kpoly[2] | at5] | + dq[5] Cpoly[3] al5] 


s=q[5] 


5 3 


y[5] 


= 0 


als] 


NSoBPII] * 2 q[5] Kpoly[3] # a[5] Kpoly[5)] 
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d 


+ 


dq[5] Cpoly[5] = 0 


3 


a[6] Kpoly[2] | al6] | 


T 


3 
gl6] Kpoly 4 alei l = dgis] Ceelv[2] | a LEI 


2 
dale] Ccpoly[4] | alo] | + dais] Gpolyi3ii gic) 


4 
CELEJ Cpoly [1] aleli dale] Cpely[5] = dle] Kpoly[5] - ule] 


5 3 
sale. Kpolvliljv 329006] Kpolyl3] 


mb[4] ddq[2] R cos(Omega t + Phi[4]) cos(q[6]) 

mb[4] ddq[2] R sin(Omega t + Phi[4]) sin(q[6]) 

mb[4] ddq[1] R sin(Omega t + Phi[4]) cos(q[6]) 

mb[4] ddq[1] R cos(Omega t + Phi[4]) sin(q[6]) 
2 2 


mb[4] Omega el R sin(q[6]) » mb[4] R ddq[6] = 0 


2 2 


mb[5] R ddq[7] + mb[5] Omega el R sin(a[7]) 


mb[5] ddq[1] R cos(Omega t + Phi[5]) sin(q[7]) 
mb[5] ddq[1] R sin(Omega t + Phi[5]) cos(q[7]) 
mb[5] ddq[2] R sin(Omega t + Phi[5]) sin(q[7]) 
mb[5] ddq[2] R cos(Omega t + Phi[5]) cos(q[7]) 


2 
dq[7] Cpoly[3] q[7] + da[7] Cpoly[4] | al7] | 


3 
se mimepoly (217 [cl TA] 


3 
qizisxpely[i4]- | c2 md ESI Kportvr20) kal] 


4 5 
aam Cpely[1] Gil mai 3 ar T KOOLI 
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3 
+ 2 ql7] Kpolvi3] * q[7] Xpolv[5] * dq[7] Cpolv[5] - 0 
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APPENDIX H - 5 BLADED POLYNOMIAL S-FUNCTION INPUT FILE 


function [I1,12,13,174,75, 16]-poly5iın 

% This m-file serves as input file for running the simulink 

% S-function simple5p.m for the 5 bladed RAH-66 Wind Tunnel Model. 
% 

EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEESEEEEEEEEEEEEEEEEEEEEEEEEEEES 

% Helo Physical and Aerodynamic parameters 
EEEEEEEEEEEEEEEEEEEEEESEEEEEESESEEEEEEEEEEEEEEEEEEEEEEEEEESES 


% Equivalent Length of rotor blade (length). 

REZO Sado % inches 

% Rotor speed (radians per/sec) 

Omega=870* (2*pi/60); 3 Run + 

% Hinge offset (length). 

el=8.5235; $ inches (9.804, 7.6898 ??) 

% Lead/lag stop position (rad). 

z=pi; 

$ Mass of rotor blades (mass). 

mb(1)20.013074; $ slugs 

mo /2)20.013074, 

mb (3)=0.013074; 

mb(4)<0.013074; 

mb(5)<0.013074; 

% Effective mass of fuselage. 

hz0.610; $ hub offset in meters 
M(1)23.504/h^2*0.06852; % convert fuselage inertia to equivalent 
M(2)27.495/h^2*0.06852; $ mass at hub & convert to slugs 

$ Blade azimuth phase angles (radians) 

Phi(1)=0; 

Phi (2)=2*pi/5; 

Phi (3) =4*p1/5; 

Ph1(4)-6:pi/5; 

Phi(5)<8"pi/5; 

% Blade Parameters 

% Spring stiffness polynomial for lead-lag (moment/radian) 


$ snubber stiffness(lb/in) * geogain(in/rad) * (ei- 
r snub) (in) 
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gg=5.1121; % geometric gain of snubber movement from lag, in/rad 
% Kpoly=[0,0,0,0,79.71*g9g+660.9141]* (le1-2.5984); % Linearized 1995 
snubber 

Kbeam-[0 00 0 660.9141]; $ original Myklestad output 
€ Ksnub=[38512436.84*gg*5,-4491019.58*gg*4,183336.70*gg*3,- 
3147:25"gg 2m 
% 19.71 *og)* (el—2. 59340. 

% 1995 snubber 

Ksnub=[336241084.22*gg*5,-38666827.45*gg*4,1547897.06*gg*3,- 
26483.84 er Du 

233. 32’cgg|*(el-2.5984), 

% 1992 snubber 
$ Kpoly=Kbeam+Ksnub; 

Kpoly=2 .5* (Kbeam+Ksnub) ; 


% Lead-lag stop spring constants (moment/radian) 


Ks (1) =0; 
Ks (2) =0; 
Ks (3) =0; 
Ks (4) =0; 
Ks (5)=0; 


$ Linear damping in lead-lag (moment/ (rad/sec) ) 
$ snubber damping(lb*s/in) * geogain(in/rad) " (el- 
r snub) (in) 


$ Cpoly=[0,0,0,0,0.59*gg]* (el-2.5984); % 
linearized 1995 snubber 
% Cpoly=[116695.68*gg*5, -14022.78*gg*4,609.0*gg*3,- 
12.05*gg^2,0.59*gg]... 
$ *(e1-2.5984); 

% 1995 snubber 

Cpoly=(-842875.01*gg"5,87576.78*gg"4,- 
20517 0] 599053225 00]*9090 2.0:47*9gg9g] 
*(e1-2.5984); 
$ 1992 snubber 
Cpolys2*Cpoly; 


$ Fuselage Parameters 


$ Linear springs in translation (force/length) 

Bey =-1216265/n22*.2248 73903 7 1227.14, % convert rotary 
spring to lateral 

E02 1-1]193.-5/h^2* 22487 3993 4912€ % spring at hub plane & 


convert to Ibf/in 


$ Linear damping in translation (force/(length/sec) ) 
$ convert rotary damping to lateral 
$ damping at hub plane & convert to lbf/in/s 
% ? c(1)225909254 0c cse T. 43/70:254-(1.15*2*p1)^2)*.06852; 
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oe ge dP 


11 


I2 


I3 


I4 


I5 


I6 


e(2)=2*3275/6 22sart (64.0573.75=- (0. 6522791) 22,7.,06852, 
c(1)=6.5820022481739737/h72,; 
C(2)=10. 2150022421* 3913 7/h 25 
cEM*2-*[0.0467,0.0542].. SSgrt CK. /M)S 


Non-linear damping in translation(force/(length/sec)^2) 


EE 
Ee 


Initial conditions 


Blade lead-lag displacement (radians) 


TITO, 
oT. 
KSAS OF 
41 = 0): 
XSS OF 


Blade lead-lag rates (radians/sec) 
EE 
xr2i=0; 
perge (OI 
xr41=0; 
xr5i=0; 


Fuselage translational displacements (length) 


xo OMIT 
zy jg pee 


Fuselage translational rates (length/sec) 


EX ISO; 
xrYi=0; 


Create input arrays: 


[mb, M]; 
[R,Omega,el,z]; 
Phi; 

le, v,.Cpoly); 
(Kpoly, Ks, K]; 


[xrxi KIN xr11,xr2i,xr 30 xr Kr 5] XI KYA KIT, XI, KG 40] X502] 


149 








function [sys, 


aga 


12 


13 


I4 


ES 


I6 


de de de de ge ge ge de ge de ge ge de dP dP ge ge de de de oe de ge de dP dP oe ge ge de ge oe oe oe 


ge 


40 


Cpoly 


Kpoly 
Ks 


oe JE de de de de ge dP dP ge dP dP oe ge 
IO 
y 
H 


APPENDIX I - 5 BLADED POLYNOMIAL S-FUNCTION 


freedom, 
and lead-lag rotor blade degrees of freedom. 


x0] ez blade5pl(t,x, 40, tElag, 1L,12 13,14 15 160) 


S-function arguments: 


time 

state vector 

input vector 

switch used by numerical integration (simulation) 
routine to access certain parts of the s-function 


S-function input parameters: 


(m1), mb (2), ,mb(3),mb(4) ,mb(5) ,M(T),M(2)] 
[R,Omega,el,z] 
[Phitl1).PHT(2),Phi(3)7Phi(4),PBhri(5)]] 


[c(1)76(2) 2w 0D) 7 v 02)5 
Cpoly (1) Cpoly(2),Cpolv(/(3) .Cpoly (2)  Gbolv t») 


[Kpoly (1) ,Kpoly (2) ,Kpoly(3) ,Kpoly (4) ,Kpoly(5), 
Ks(1),Ks(2),Ks(3),Ks(4),Ks(5), 
EC) K(2)] 


Dera or Ya xr li .Xr 21, XE31 xr Xr HI, 
Sa KI SA Kol, ead, MOL! 


S-function to represent dynamics of 5 bladed coupled rotor- 
fuselage model which considers only inplane degrees of 


i.e., x and y translational fuselage degrees of freedom 


Explaination of variables: 


€— mm ep zm zm mm mm mm zm zm mm zm zm mm zm vm vm vm mm mm vm zm mmm pn — 


mass of blade 

effective mass of fuselage 

distance from lead-lag hinge to blade center of mass 
blade hinge offset 

rotor speed 

angle at which blade hits stops 

blade phase angle w.r.t. azimuth postion 

fuselage linear damping 

fuselage hydraulic damping 

blade damping polynomial coefficients 

effective stiffness of fuselage (landing gear stiffness) 
blade elastic spring constant polynomial coefficients 
blade stop effective spring constant 
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% xr i Se initial rate 

% x__i -> initial displacement 

% 
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEESEEEEEES 
% Define input parameters 
ESEEESEEEEEEESESESEEEESEESEEEEEEEEEESEEEESEEEEEEESEEFEEEEEEEEEESEEEES 


mb<I1(1:5); 

MTI) 

R=I2 (1 Omega=12 (2) e1l=12(3),2=12 (4) 

Phi=13; 

C=14 (iis 2 

= A TE 

Cpoly=14 (5:9); 

Kpoly=15 (1:5); 

Ks=15(6:10); 

K=15 (11:12); 

KEX ISMON IDE XTY IS NOZE 
xr412$603)5Xr27-T6(4);xr31s16(5);xr4isTO(6)52x5522 60 0M 
2x -ISUB) yi-T6(9). 

SI 76(10), 221-176 (11) +9 1=16(12)6 241 =1:6 (013) ;x5%=16 (1475 


Ss SR 555555555233 5555 5355 35555 
% S-function flag conditionals 
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 
if flag == 


syvs=[14,0,14,7,0,0]; 
KOS TY! xr Kee ss srapocr51 8%X3,%XY),x11,X21,%X31,x41, 51]. 
elseif abs(flag) == 1 


t2 = mb(1)*R; 

t3 = Omega*t; 

t4 = t3+Phi (1); 
t5 s cos(t4); 

t6 = sin(xt10)); 
te =| to te: 

t9 - sin(t4); 

Emo cost(x(10)); 
peki t9*t10; 
i= bt yee Ls 
t14 = mb(2)*R; 
fi5 = E3+PR1 (2): 
Giles ="cos( tls); 


EI: SIN (Xa): 

cle 2 tlo t17: 

E20 z SIN (t15): 

pol: cos x( IM), 

E22. - :20*E2]5 

t24 = -t14*t18-t14*t22; 


belo ee Mota) Rs 
E26 = t3tPBh2(03)5 
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227. = cos (625) 

t28 = sin(x (129); 

t29 = t27°t28; 

t31 = sinlt26); 

E32 = COS (NG 

E335 t3l565225 

t35.2 -t25*t29-t25*E335 

t36 - mb(4)*R; 

EST. - t3rPnhi(4); 

t38 - cos(t37); 

BORS ESMASE 

t40 = t38*t39; 

ET s Ssin(u37); 

CAST cos (x(13)): 

t44 = t42*t43; 

t46 = -t36*t40-t36*t44; 

t47 = mb(5)*R; 

£48 = t3+Ph1 (5); 

t49 - cos(t48); 

t50e— sin (x (14 e 

E51 = t49*t50; 

too) Sin (tac); 

t54 = cos(x(14)); 

foo — t53*t54; 

C57 2 -t47*t51-t47*56555; 

t62 - Omega^2; 

C63 EE E52: 

t66 # C20 C17: 

t67 z x(4)72; 

E70 = C25*ts i; 

t72 = Omega*t28*x(5); 

t74 = t42*t39; 

t75 2 x(6)^2; 

t78 = t42*t62; 

E81 = t38*t62; 

t84 = t38*t43; 

te s tb356550; 

wee x2; 

E91. t36*t38; 

t93 = Omega*t43*x (6); 

t95 = t360*t42; 

t97 = Omega*t39*x(6); 

t99 = t49*t54; 

t102 2 E53*5852; 

E105 = t16*e62; 

ELOS =" £2765; 

t110 = Omega*t10*x(3); 

tama. = t9*t5; 

E113. x(3)52; 

EE 

EE 

t 1220 = =v (1) *x (1) *abs4x (27) =K(1) *x (8) +e25*e63° ES 251 toc eos > 
2*t70*t72-t36*t74*t75-t36*t78*t359+t36*t81*t43+t36*t84*te75- 
t47*t87*t88+2*t91*t93-2*t95*t97+t47*t99*t88- 
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4456 102:t50 +614 105*t21+2*t108*t110-— 
2411120 CS +6477.11116.*t54+mb(4)“Ee38*t120:; 

t131 = €47*t49; 

t133 = Omega*t54*x(7); 

t 135. = r20*t62; 

E1382=5€47*t653; 

t140 = Omega*t50*x(7); 

tl42 - tl6*E21:; 

t145 = t14*t16; 


t147 = Omega*t21*x(4); 
t149 = t14*t20; 
t151 = Omega*t17*x (4); 


E53 t0 t62: 

ti5o6Um t5*t62;: 

EL ti1*t28; 

E 1605=.x (51 2; 

plošče ES ELO; 

tl66.— t2*t9; 

t168 = Omega*t6*x(3); 

E 171 2 46234632, 

674 = t31*t62:; 

El GN kA]: 

E179 s Omega*t32*x(5); 

t181 - 
55149 12 0 m5(1*tE5*t120+mB (2) *=16*t 120585) *t2142120+2*t131*t 1535 
146 185*e17—-2*E138*t140+t14*t142*t67t2*t145*t147+u(1)-2*t149*t151= 
t2*Ll1l53*t6t*t2*t156*t10-t25*t159*t160*t2*t163*t113-2*t166*Ct168- 
ci1)*x001)*6525*E171*tE160-t25*t174*t2842*t177*t179; 

e186 - -t2*tll2*t2*tl153; 

t189 = t14*t142-t14*t66; 


E1072 -E255 6159 +E25*E 171; 
E155 =5636*E74+636*E84; 
t198 = -t47*t87+t47*t99; 


EE - 
EE EE EE EE EE 
GE T2 
e oe EE 7k95 TENG 
EE EE 

t265 z 
t47g*t5lsLOSBTEd7*tllo*ct5O0«*2*tl45*tl151ttE2*C7* E WIE 2*C 1 38B*E T3353 E47 ^ E555 E SS 
+261 51541102 e166*C110+t2*E11*E113+t47*t1024t54+t14e105*E17+814*2 2 
Lto7zxmbCtOSETZOFEIA*EI135*t21-«mbD(4)*t42*t120*mD(5)*E53*t120*mb(3) *t 34 
t 12052. E 108168=K(2) *x (9) +C2*EC153*E10; 


01267 2 R^2; 
5269 2 x1(110)^2; 
EDIT 205 


t274 = abs(x(10)); 

E275 = 121472; 

E2765 ez et 056214: 

E287 2 626972: 

t299 = el*R; 

1802 = =t2698Koolv (4) *E271/2=x (10) *Kpoly (2) *t276- 
x lO KpBo1lv. (4) 2224 =x 13) "Cpoly (2) E276=x (3) *Cpolvyv (4) *E274- 
Rio) Epa UY ST 2269=x (3) *Cpoly (8) t267+4u (3) = 


154 


3.E0/2.E0*t269*Kpoly (2) *t275*t271-x(10) *Kpoly (5) -x(3) *Cpoly (5) - 
3*t287*x(10) *Kpoly (1) -2*t269*x(10) *Kpoly (3) -mb(1) *t62*t299*t6; 
ESU5. —c TS 


E306 - L5UB 2; 
t318 - abs(x(11)); 
65 194 = 1531802 
t320 = 10; 

E324 2 E319*t318; 


t336 - u(4)-x(4) *Cpoly(1)*t306-3*t306*x(11) *Kpoly(1)- 
Q*t305*x(11)*Kpoly(3)-x(11)*Kpoly(5)-x(4)*Cpoly(5)-mb(2)*t62*8299*t17- 
3.50/2.E0*t305*Kpoly(2)*t319*t320-x(4) *Cpoly(2)*t924- 
x(11)*Kpoly(2)*t324-x(11)*Kpoly(4)*t318-t305*Kpoly(4) *t32072- 
x(4)*Cpolv(3)*6E305-x(4)*Cpoly(4)*t318; 

Page =o mee) 25 

E310 2 £3392; 

£545) = abstx(12))- 

£348 = t345^2; 

t349 = t348*t345; 

t356 = 0; 

t370 = u(5)-x(5) *Cpoly (1) *t340-x(5) *Cpoly (3) *t339- 
x(5) *Cpoly (4) *t345-x(12) *Kpoly (2) *t349-x (12) *Kpoly (4) “t345- 
(5) *Cpoly(2)*t349-=t339*Kpoly (4) *t356/2-2*t339*x(12) "Kpoly(3)= 
3.5E0/2.E0*t339*Kpoly(2)*t348*t356-mb(3)*t62*t299*t28=Xx(5) *Cpoly (5) -= 
x(12)*Kpoly(5)-3*t340*x(12)*Kpoly(1); 


EO Ex (l3) 2- 
t374 = abs (x(13)); 
S750 = 0574 2° 
t376 2.0: 

ESSI < E5/202; 


E390 =7E3 7576374: 

E404 - ut6)-3.E07/2:R0*tE372*Kpolv(2)*E375*t376-23: (13) *REOS 5) a 
paU6)*Cpoly(5)-3*t381*x(13)*Kpoly(1)-2*t372*x(13)*Kpolw(3)- 
mb(4)*t62*t299*t39-x(6) *Cpoly(2)*t390-x(6)*Cpoly(4)*t374- 
x(13)*Kpoly(4)*t374-t372*Kpoly(4)*t376/2-x(13)*Rpoly(2)*t390- 
2461 CBolv (3) *E372=x(6)"Cpoly (1) *E381; 

t406 e x(14)^2; 

t408 = abs(x(14)); 

t409 = t408°2; 

tato = 0; 

t421 = t409*t408; 

t431 = t406°2; 

t438 = u(7)-3.E0/2.E0*t406*Kpoly (2) *t409*t410-x (7) *Cpoly (4) *t408- 
mb(5)*t62*t299*t50-t406*Kpoly(4)*t410/2-x(Y14)*Kpoly(2)*E421—- 
X(14) "Kpoly(4) "t408-x(7) "Cpoly(2) "t421-x(7) "Cpoly(3) "t406- 
x(7)*Cpoly(5)-x(14) *Kpoly(5)-3*t431*x(14)*Kpoly(1)- 

PAREADO Xx (14) *kKpoly(3)=x(7) *Cpoly (1) *t431; 

B(1,1) = mb(4)+mb(5)+mb(1)+mb(2)+mb(3)+M(1); 

BUM 2) = O- 

EN 200 = 137. 

Bli;4) = C24: 

Bili 5m et 35: 

S(O) S ede} 

Bi a) C51: 

B(1,8) = t122+t181; 
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B( 27.1) 
B(2,2) 
Bl) 
B2 4 
B(2,5) 
B(2,6) 
E27) 
B(2,8) 
B(3,1) 
B(3,2) 
B(3,3) 
B(3,4) 
B3 5) 
B(3,6) 
B(3,7) 
B(3,8) 
B(4 1) 
B(4,2) 
B(4,3) 
B(4,4) 
B(4,5) 
B(4,6) 
B (477) 
B(4,8) 
B (571) 
B(5,2) 
B5, 3) 
B(5,4) 
B(5,5) 
B(5,6) 
B(5,7) 
B(5,8) 
B(6,1) 
B(6,2) 
B(6,3) 
B(6,4) 
B(6,5) 
B(6,6) 
B(6,7) 
B(6,8) 
B(7,1) 
BIGZ) 
B(7,3) 
B(7,4) 
B(7,5) 
B(7,6) 
B(7,7) 
B(7,8) 


OF 


= mb(1)+mb (2) +mb (3) +mb (4) +mb (5) +M (2); 
SES 


t189; 

E192> 

E195. 

t198, 
23046265, 
E136 

t186; 
m5(1)*t267; 
0; 


“2267; 


= mb(4) *t267; 


H 


O; 
mb (5) 
t438; 


E SZG 


% Calculate derivatives 


[m,n]ssize(B); 
ATSB(:,l:m-l); 
EB: rn), 
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sys=zeros (1.2 :m)E 
sys(1:7)=A11£1; 
sys (8:14)=x(1:7); 
% Output states 
elseif abs(flag) == 3 
Sys sx: 
else 
sys = []; 


end 
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